rm(list = ls())

#-----Setup------
suppressMessages({
  library(tidyverse) # for easy data cleaning
  library(lubridate) # for working w/ dates
  library(estimatr)
  library(coefplot)
  library(ggthemes)
  library(lemon)
  library(ggforce)
  library(ggrepel)
  library(scales)
  library(grid)
  library(cowplot)
  library(haven)
  library(ggbeeswarm)
})


## Helper functions:
add_parens <- function(x, digits = 2) {
  x <- as.numeric(x)
  return(paste0("(", sprintf(paste0("%.", digits, "f"), x), ")"))
}

format_num <- function(x, digits = 2) {
  x <- as.numeric(x)
  return(paste0(sprintf(paste0("%.", digits, "f"), x)))
}

make_entry <- function(est, se, p, digits = 2) {
  entry <- paste0(format_num(est, digits = digits),
                  " ",
                  add_parens(se, digits = digits))
  entry[p < 0.05] <- paste0(entry[p < 0.05], "*")
  return(entry)
}

table_entry <- function(est, se, digits = 2) {
  entry <- paste0(format_num(est, digits = digits),
                  " ",
                  add_parens(se, digits = digits))
  return(entry)
}


### ------ Import data ------
dat <- read_dta("lucid_oct29_clean.dta")

## Set plot parameters
bpr_colors <- c("#1F3A93", "#7C2C55", "#D91E18")

fixed <- 2
spacing <- .285

## Pre-treatmenet covariates for regression adjustment
covariates <-
  c(
    "dem_age_cat",
    "dem_sex_cat",
    "dem_region_cat",
    "dem_hhi_cat",
    "dem_emp_cat",
    "dem_race_cat",
    "dem_edu_cat",
    "dem_pid_7n",
    "pol_ideo_7n",
    "pol_attent_n"
  )

## Recodes for easier interpretation of estimates
dat <-
  dat %>%
  mutate(
    Z = factor(Z, levels = c(1:3), labels = c("Reform", "Defund", "Abolish")),
    ## Make these binary for oppose rather than support
    othr_act_violent_n   = 1 - othr_act_violent_n,
    othr_act_property_n  = 1 - othr_act_property_n,
    othr_act_vice_n      = 1 - othr_act_vice_n,
    othr_act_order_n     = 1 - othr_act_order_n,
    othr_act_prevent_n   = 1 - othr_act_prevent_n,
    othr_act_traffic_n   = 1 - othr_act_traffic_n,
    othr_act_social_n    = 1 - othr_act_social_n,
    othr_act_engage_n    = 1 - othr_act_engage_n,
    ## Recode to run from increase to decrease
    othr_change_numcops_n = as.numeric(6 - othr_change_numcops_n),
    othr_spend_police_n = as.numeric(6 - othr_spend_police_n),
    othr_spend_edu_n     = as.numeric(6 - othr_spend_edu_n),
    othr_spend_health_n  = as.numeric(6 - othr_spend_health_n),
    othr_spend_housing_n = as.numeric(6 - othr_spend_housing_n),
    othr_spend_social_n  = as.numeric(6 - othr_spend_social_n),
    ## Reverse code and create binary DVs to estimate proportions
    cuts_safe_n = as.numeric(6 -  cuts_safe_n),
    cuts_safe_binary = as.numeric(cuts_safe_n > 3),
    cuts_crime_binary = as.numeric(cuts_crime_n > 3),
  ) 

### ------ Fig. 1 Manuscript -------

### Within-subjects estimates for labels v. substance.
gg_u_df <-
  dat %>%
  select(supp_reform_n,
         supp_defund_n,
         supp_abolish_n) %>%
  pivot_longer(supp_reform_n:supp_abolish_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(DV) %>%
  do(tidy(lm_robust(Y ~ 1, data = .))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    DV = case_when(
      DV == "supp_reform_n" ~  "Reforming the police",
      DV == "supp_defund_n" ~  "Defunding the police",
      DV == "supp_abolish_n" ~ "Abolishing the police"
    ),
    type = "Slogan"
  ) %>%
  bind_rows(
    dat %>%
      select(policy_reform_n,
             policy_defund_n,
             policy_abolish_n) %>%
      pivot_longer(
        policy_reform_n:policy_abolish_n,
        names_to = "DV",
        values_to = "Y"
      ) %>%
      group_by(DV) %>%
      do(tidy(lm_robust(Y ~ 1, data = .))) %>%
      filter(term == "(Intercept)") %>%
      mutate(
        DV = case_when(
          DV == "policy_reform_n" ~  "Reforming the police",
          DV == "policy_defund_n" ~  "Defunding the police",
          DV == "policy_abolish_n" ~ "Abolishing the police"
        ),
        type = "Substance"
      )
  )


gg_w_df <-
  dat %>%
  select(supp_reform_n,
         supp_defund_n,
         supp_abolish_n,
         weights) %>%
  pivot_longer(supp_reform_n:supp_abolish_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(DV) %>%
  do(tidy(lm_robust(Y ~ 1, weights = weights,  data = .))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    DV = case_when(
      DV == "supp_reform_n" ~  "Reforming the police",
      DV == "supp_defund_n" ~  "Defunding the police",
      DV == "supp_abolish_n" ~ "Abolishing the police"
    ),
    type = "Slogan"
  ) %>%
  bind_rows(
    dat %>%
      select(policy_reform_n,
             policy_defund_n,
             policy_abolish_n,
             weights) %>%
      pivot_longer(
        policy_reform_n:policy_abolish_n,
        names_to = "DV",
        values_to = "Y"
      ) %>%
      group_by(DV) %>%
      do(tidy(lm_robust(
        Y ~ 1, data = ., weights = weights
      ))) %>%
      filter(term == "(Intercept)") %>%
      mutate(
        DV = case_when(
          DV == "policy_reform_n" ~  "Reforming the police",
          DV == "policy_defund_n" ~  "Defunding the police",
          DV == "policy_abolish_n" ~ "Abolishing the police"
        ),
        type = "Substance"
      )
  )

gg_df <-
  bind_rows(
    gg_u_df %>% mutate(estimator = "Unweighted"),
    gg_w_df %>% mutate(estimator = "Weighted")
  ) %>%
  mutate(type = factor(type, levels = c("Substance", "Slogan")))

g <-
  gg_df %>%
  filter(estimator == "Weighted") %>%
  ggplot(.,
         aes(
           x = estimate,
           y = type,
           color = estimator,
           fill = estimator,
           group = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.5),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.5),
    size = 3,
    pch = 21,
    color = "white"
  ) +
  scale_fill_manual("", values = c("black")) +
  scale_color_manual("", values = c("black")) +
  facet_col( ~ DV) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    #strip.text=element_text(hjust = 0.55),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Respondents' average level of support (1 = None at all,...,5 = A great deal)",
       y = "")

g <- g + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g

g %>%
  ggsave(
    filename = "Fig1.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * (nrow(g$data) * .75)
  )

### ------ Fig. B1 Appendix ----

## Between-subjects estimates for labels v. substance.
gg_u_df <-
  dat %>%
  filter(Z_label_first == 1) %>%
  select(supp_reform_n,
         supp_defund_n,
         supp_abolish_n) %>%
  pivot_longer(supp_reform_n:supp_abolish_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(DV) %>%
  do(tidy(lm_robust(Y ~ 1, data = .))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    DV = case_when(
      DV == "supp_reform_n" ~  "Reforming the police",
      DV == "supp_defund_n" ~  "Defunding the police",
      DV == "supp_abolish_n" ~ "Abolishing the police"
    ),
    type = "Slogan"
  ) %>%
  bind_rows(
    dat %>%
      filter(Z_label_first == 0) %>%
      select(policy_reform_n,
             policy_defund_n,
             policy_abolish_n) %>%
      pivot_longer(
        policy_reform_n:policy_abolish_n,
        names_to = "DV",
        values_to = "Y"
      ) %>%
      group_by(DV) %>%
      do(tidy(lm_robust(Y ~ 1, data = .))) %>%
      filter(term == "(Intercept)") %>%
      mutate(
        DV = case_when(
          DV == "policy_reform_n" ~  "Reforming the police",
          DV == "policy_defund_n" ~  "Defunding the police",
          DV == "policy_abolish_n" ~ "Abolishing the police"
        ),
        type = "Substance"
      )
  )


gg_w_df <-
  dat %>%
  filter(Z_label_first == 1) %>%
  select(supp_reform_n,
         supp_defund_n,
         supp_abolish_n,
         weights) %>%
  pivot_longer(supp_reform_n:supp_abolish_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(DV) %>%
  do(tidy(lm_robust(Y ~ 1, weights = weights,  data = .))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    DV = case_when(
      DV == "supp_reform_n" ~  "Reforming the police",
      DV == "supp_defund_n" ~  "Defunding the police",
      DV == "supp_abolish_n" ~ "Abolishing the police"
    ),
    type = "Slogan"
  ) %>%
  bind_rows(
    dat %>%
      filter(Z_label_first == 0) %>%
      select(policy_reform_n,
             policy_defund_n,
             policy_abolish_n,
             weights) %>%
      pivot_longer(
        policy_reform_n:policy_abolish_n,
        names_to = "DV",
        values_to = "Y"
      ) %>%
      group_by(DV) %>%
      do(tidy(lm_robust(
        Y ~ 1, data = ., weights = weights
      ))) %>%
      filter(term == "(Intercept)") %>%
      mutate(
        DV = case_when(
          DV == "policy_reform_n" ~  "Reforming the police",
          DV == "policy_defund_n" ~  "Defunding the police",
          DV == "policy_abolish_n" ~ "Abolishing the police"
        ),
        type = "Substance"
      )
  )


gg_robust <-
  bind_rows(
    gg_u_df %>% mutate(estimator = "Unweighted"),
    gg_w_df %>% mutate(estimator = "Weighted")
  ) %>%
  mutate(type = factor(type, levels = c("Substance", "Slogan")),
         estimand = "Between-subjects")

gg_df <-
  bind_rows(gg_df %>% mutate(estimand = "Within-subjects"),
            gg_robust)

g1 <-
  gg_df %>%
  filter(estimand == "Within-subjects") %>%
  ggplot(.,
         aes(
           x = estimate,
           y = type,
           fill = estimator,
           group = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 2
  ) +
  geom_point(position = position_dodgev(height = 0.7),
             size = 3,
             pch = 21) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  facet_col( ~ DV) +
  scale_x_continuous(limits = c(1, 5)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "",
       y = "")

g2 <-
  gg_df %>%
  filter(estimand != "Within-subjects") %>%
  ggplot(.,
         aes(
           x = estimate,
           y = type,
           fill = estimator,
           group = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 2
  ) +
  geom_point(position = position_dodgev(height = 0.7),
             size = 3,
             pch = 21) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  facet_col( ~ DV) +
  scale_x_continuous(limits = c(1, 5)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Respondents' average level of support\n (1 = None at all,...,5 = A great deal)",
       y = "")

g1 <- g1 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g2 <- g2 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))

g <-
  plot_grid(
    g1,
    g2,
    ncol = 1,
    labels = "auto",
    label_fontface = "plain",
    label_fontfamily = "sans",
    label_size = 12
  )

g

g %>%
  ggsave(filename = "FigB1.pdf",
         .,
         width = 6.5,
         height = 7)

### ------ Table B1.1 Appendix -------------

tab_top <-
  g1$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -term, -.group, -outcome) %>%
  arrange(estimator, DV, desc(type)) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, type, entry) %>%
  pivot_wider(names_from = type:estimator,
              values_from = entry)

tab_btm <-
  g2$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -term, -.group, -outcome) %>%
  arrange(estimator, DV, desc(type)) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, type, entry) %>%
  pivot_wider(names_from = type:estimator,
              values_from = entry)

bind_rows(tab_top, tab_btm) %>%
  select(
    DV,
    `Slogan_Unweighted`,
    `Slogan_Weighted`,
    `Substance_Unweighted`,
    `Substance_Weighted`
  ) %>%
  write.table(
    .,
    file = "TabB1.1.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

### ------ Table B1.2 Appendix ----------

### Within-subjects estimates for labels v. substance.

gg_df <-
  dat %>%
  select(
    weights,
    covariates,
    admin_ResponseId,
    supp_reform_n,
    supp_defund_n,
    supp_abolish_n,
    policy_reform_n,
    policy_defund_n,
    policy_abolish_n
  ) %>%
  pivot_longer(supp_reform_n:policy_abolish_n,
               names_to = "Z",
               values_to = "Y") %>%
  mutate(
    Z_movement = case_when(
      str_detect(Z, "reform") ~ "Reform",
      str_detect(Z, "defund") ~ "Defund",
      str_detect(Z, "abolish") ~ "Abolish"
    ),
    Z_type = case_when(
      str_detect(Z, "supp") ~ "Slogan",
      str_detect(Z, "policy") ~ "Policy"
    ),
    Z_type = factor(Z_type, levels = c("Slogan", "Policy"))
  )

## Difference in means.
gg_dim_u <-
  gg_df %>%
  group_by(Z_movement) %>%
  do(tidy(
    lm_robust(
      Y ~ Z_type,
      clusters = admin_ResponseId,
      se_type = "stata",
      data = .
    )
  )) %>%
  mutate(estimator = "Diff-in-Means",
         weighted = FALSE)

gg_dim_w <-
  gg_df %>%
  group_by(Z_movement) %>%
  do(tidy(
    lm_robust(
      Y ~ Z_type,
      clusters = admin_ResponseId,
      se_type = "stata",
      weights = weights,
      data = .
    )
  )) %>%
  mutate(estimator = "Diff-in-Means",
         weighted = TRUE)

## Covariate-adjusted
gg_adj_u <-
  gg_df %>%
  group_by(Z_movement) %>%
  do(tidy(
    lm_lin(
      Y ~ Z_type,
      covariates = as.formula(paste("~", paste(
        covariates, collapse = "+"
      ))),
      clusters = admin_ResponseId,
      se_type = "stata",
      data = .
    )
  )) %>%
  mutate(estimator = "Covariate-adjusted",
         weighted = FALSE)

gg_adj_w <-
  gg_df %>%
  group_by(Z_movement) %>%
  do(tidy(
    lm_lin(
      Y ~ Z_type,
      covariates = as.formula(paste("~", paste(
        covariates, collapse = "+"
      ))),
      clusters = admin_ResponseId,
      se_type = "stata",
      weights = weights,
      data = .
    )
  )) %>%
  mutate(estimator = "Covariate-adjusted",
         weighted = TRUE)

## Combine estimates
est_df <-
  bind_rows(gg_dim_u, gg_dim_w, gg_adj_u, gg_adj_w) %>%
  select(-outcome, -df) %>%
  filter(term %in% c("Z_typePolicy", "(Intercept)")) %>%
  mutate(term = case_when(
    term == "(Intercept)" ~ "Slogan Mean",
    term == "Z_typePolicy" ~ "Policy Effect"
  ))

## Generate table for word
ols_table <-
  est_df %>%
  filter(term != "Slogan Mean") %>%
  mutate(entry = make_entry(
    est = estimate,
    se = std.error,
    p = p.value,
    digits = 2
  )) %>%
  dplyr::select(estimator, weighted, entry) %>%
  pivot_wider(names_from = estimator:weighted,
              values_from = c(entry))


write.table(
  ols_table,
  file = "TabB1.2.txt",
  sep = ",",
  quote = FALSE,
  row.names = F,
  col.names = FALSE
)


### ------ Estimate group means -------
gg_u_df <-
  dat %>%
  select(
    Z,
    othr_supp_eliminate_n,
    othr_supp_replace_n,
    othr_spend_police_n,
    othr_change_numcops_n,
    othr_supp_destroy_n,
    othr_supp_attack_n,
    othr_supp_verbal_n,
    othr_supp_choke_n,
    othr_supp_knock_n,
    othr_supp_military_n,
    othr_supp_identify_n,
    othr_supp_crb_n,
    othr_supp_bwc_n,
    othr_supp_live_n,
    othr_supp_protest_brutal_n,
    othr_supp_protest_blm_n,
    othr_supp_protest_defund_n,
    othr_act_violent_n,
    othr_act_property_n,
    othr_act_vice_n,
    othr_act_order_n,
    othr_act_prevent_n,
    othr_act_traffic_n,
    othr_act_social_n,
    othr_act_engage_n,
    othr_spend_edu_n,
    othr_spend_health_n,
    othr_spend_housing_n,
    othr_spend_social_n
  ) %>%
  pivot_longer(
    othr_supp_eliminate_n:othr_spend_social_n,
    names_to = "DV",
    values_to = "Y"
  ) %>%
  group_by(Z, DV) %>%
  do(tidy(lm_robust(Y ~ 1, data = .))) %>%
  filter(term == "(Intercept)") %>%
  mutate(estimator = "Unweighted")

gg_w_df <-
  dat %>%
  select(
    Z,
    weights,
    othr_supp_eliminate_n,
    othr_supp_replace_n,
    othr_spend_police_n,
    othr_change_numcops_n,
    othr_supp_destroy_n,
    othr_supp_attack_n,
    othr_supp_verbal_n,
    othr_supp_choke_n,
    othr_supp_knock_n,
    othr_supp_military_n,
    othr_supp_identify_n,
    othr_supp_crb_n,
    othr_supp_bwc_n,
    othr_supp_live_n,
    othr_supp_protest_brutal_n,
    othr_supp_protest_blm_n,
    othr_supp_protest_defund_n,
    othr_act_violent_n,
    othr_act_property_n,
    othr_act_vice_n,
    othr_act_order_n,
    othr_act_prevent_n,
    othr_act_traffic_n,
    othr_act_social_n,
    othr_act_engage_n,
    othr_spend_edu_n,
    othr_spend_health_n,
    othr_spend_housing_n,
    othr_spend_social_n
  ) %>%
  pivot_longer(
    othr_supp_eliminate_n:othr_spend_social_n,
    names_to = "DV",
    values_to = "Y"
  ) %>%
  group_by(Z, DV) %>%
  do(tidy(lm_robust(Y ~ 1, data = ., weights = weights))) %>%
  filter(term == "(Intercept)") %>%
  mutate(estimator = "Weighted")


gg_df  <-
  bind_rows(gg_u_df, gg_w_df) %>%
  mutate(
    DV = factor(
      DV,
      levels = c(
        "othr_supp_destroy_n",
        "othr_supp_attack_n",
        "othr_supp_eliminate_n",
        "othr_supp_replace_n",
        "othr_spend_police_n",
        "othr_change_numcops_n",
        "othr_supp_military_n",
        "othr_supp_identify_n",
        "othr_supp_crb_n",
        "othr_supp_live_n",
        "othr_supp_bwc_n",
        "othr_supp_choke_n",
        "othr_supp_knock_n",
        "othr_supp_verbal_n",
        "othr_supp_protest_brutal_n",
        "othr_supp_protest_blm_n",
        "othr_supp_protest_defund_n",
        "othr_act_violent_n",
        "othr_act_property_n",
        "othr_act_vice_n",
        "othr_act_order_n",
        "othr_act_prevent_n",
        "othr_act_traffic_n",
        "othr_act_social_n",
        "othr_act_engage_n",
        "othr_spend_edu_n",
        "othr_spend_health_n",
        "othr_spend_housing_n",
        "othr_spend_social_n"
      ),
      labels = c(
        "Destruction of property",
        "Attacks on police officers",
        "Eliminating police departments",
        "Replacing police with social services",
        "Spending on police services",
        "Number of police on the street",
        "Reduce use of military equipment",
        "Make firing police easier",
        "Implement Civilian Review Boards",
        "Require police to live where they work",
        "Require Body Worn Cameras",
        "Ban use of chokeholds",
        "Ban use of no-knock warrants",
        "Require warning before shooting",
        "Protest against Police Brutality",
        "Protests by Black Lives Matter",
        "Protests to Defund the Police",
        "Responding to violent crimes",
        "Responding to property crimes",
        "Responding to vice crimes",
        "Conducting vehicle stops",
        "Crime prevention",
        "Maintaining civil order",
        "Providing social services",
        "Community engagement",
        "Spending on education",
        "Spending on healthcare",
        "Spending on housing",
        "Spending on social services"
      )
    ),
    Z = case_when(
      Z == "Abolish" ~ "Abolishing",
      Z == "Defund" ~ "Defunding",
      Z == "Reform" ~ "Reforming"
    ),
    Z = factor(Z, levels = rev(
      c("Abolishing", "Defunding", "Reforming")
    ))
  )

### ------ Fig. 2 Manuscript ----------
g1 <-
  gg_df %>%
  filter(
    DV %in% c("Destruction of property",
              "Attacks on police officers"),
    estimator == "Weighted"
  ) %>%
  ggplot(.,
         aes(x = estimate,
             y = Z)) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.7),
    size = 3,
    pch = 21,
    color = "white",
    fill = "black"
  ) +
  facet_row( ~ DV) +
  scale_x_continuous(labels = percent,
                     expand = c(0, 0),
                     limits = c(0.1, 0.9)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Percentage of respondents that believe most supporters favor this",
       y = "")


## Read in data from supplementary experiment
dat2 <-
  read_dta("lucid_feb18_clean.dta") %>%
  mutate(Z = factor(
    Z,
    levels = c(1:3),
    labels = c("Reform", "Defund",
               "Abolish")
  ))

gg_u_df2 <-
  dat2 %>%
  select(Z,
         weights,
         starts_with("Y_")) %>%
  pivot_longer(Y_black_support_n:Y_voter_support_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(Z, DV) %>%
  do(tidy(lm_robust(Y ~ 1, data = .))) %>%
  mutate(
    reference = case_when(
      DV == "Y_black_support_n" ~ 13,
      DV == "Y_female_support_n" ~  50,
      DV == "Y_biden_support_n" ~   51,
      DV == "Y_voter_support_n" ~   67
    ),
    DV = case_when(
      DV == "Y_black_support_n" ~ "Proportion Black (reference: 13% of U.S. population)",
      DV == "Y_female_support_n" ~  "Proportion Female (reference: 50% of U.S. population)",
      DV == "Y_biden_support_n" ~   "Proportion voted for Biden (reference: 51% of voting population)",
      DV == "Y_voter_support_n" ~   "Proportion voted in 2020 (reference: 67% of eligible population)"
    ),
    Z = case_when(
      Z == "Abolish" ~ "Abolishing",
      Z == "Defund" ~ "Defunding",
      Z == "Reform" ~ "Reforming"
    ),
    Z = factor(Z, levels = rev(
      c("Abolishing", "Defunding", "Reforming")
    ))
  )

gg_w_df2 <-
  dat2 %>%
  select(Z,
         weights,
         starts_with("Y_")) %>%
  pivot_longer(Y_black_support_n:Y_voter_support_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(Z, DV) %>%
  do(tidy(lm_robust(Y ~ 1, data = ., weights = weights))) %>%
  mutate(
    reference = case_when(
      DV == "Y_black_support_n" ~ 13,
      DV == "Y_female_support_n" ~  50,
      DV == "Y_biden_support_n" ~   51,
      DV == "Y_voter_support_n" ~   67
    ),
    DV = case_when(
      DV == "Y_black_support_n" ~ "Proportion Black (reference: 13% of U.S. population)",
      DV == "Y_female_support_n" ~  "Proportion Female (reference: 50% of U.S. population)",
      DV == "Y_biden_support_n" ~   "Proportion voted for Biden (reference: 51% of voting population)",
      DV == "Y_voter_support_n" ~   "Proportion voted in 2020 (reference: 67% of eligible population)"
    ),
    Z = case_when(
      Z == "Abolish" ~ "Abolishing",
      Z == "Defund" ~ "Defunding",
      Z == "Reform" ~ "Reforming"
    ),
    Z = factor(Z, levels = rev(
      c("Abolishing", "Defunding", "Reforming")
    ))
  )

gg_df2 <-
  bind_rows(
    gg_u_df2 %>% mutate(estimator = "Unweighted"),
    gg_w_df2 %>% mutate(estimator = "Weighted")
  )

## underlying distribution
raw_df <-
  dat2 %>%
  select(Z,
         starts_with("Y_")) %>%
  pivot_longer(Y_black_support_n:Y_voter_support_n,
               names_to = "DV",
               values_to = "Y") %>%
  mutate(
    reference = case_when(
      DV == "Y_black_support_n" ~ 13,
      DV == "Y_female_support_n" ~  50,
      DV == "Y_biden_support_n" ~   51,
      DV == "Y_voter_support_n" ~   67
    ),
    DV = case_when(
      DV == "Y_black_support_n" ~ "Proportion Black (reference: 13% of U.S. population)",
      DV == "Y_female_support_n" ~  "Proportion Female (reference: 50% of U.S. population)",
      DV == "Y_biden_support_n" ~   "Proportion voted for Biden (reference: 51% of voting population)",
      DV == "Y_voter_support_n" ~   "Proportion voted in 2020 (reference: 67% of eligible population)"
    ),
    Z = case_when(
      Z == "Abolish" ~ "Abolishing",
      Z == "Defund" ~ "Defunding",
      Z == "Reform" ~ "Reforming"
    ),
    Z = factor(Z, levels = rev(
      c("Abolishing", "Defunding", "Reforming")
    ))
  )

g2 <-
  gg_df2 %>%
  filter(estimator == "Weighted",
         DV == "Proportion Black (reference: 13% of U.S. population)") %>%
  ggplot(.,
         aes(
           x = estimate / 100,
           y = Z,
           fill = estimator,
           group = estimator
         )) +
  geom_quasirandom(
    data = raw_df %>%
      filter(DV == "Proportion Black (reference: 13% of U.S. population)"),
    aes(x = Y / 100, y = Z),
    color = "darkgrey",
    alpha = 0.4,
    inherit.aes = FALSE,
    groupOnX = FALSE,
    varwidth = TRUE
  ) +
  geom_vline(aes(xintercept = reference / 100),
             linetype = 1,
             size = 0.5) +
  geom_errorbarh(
    aes(xmin = conf.low / 100, xmax = conf.high / 100),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.7),
    size = 3,
    pch = 21,
    color = "white",
    fill = "black"
  ) +
  facet_col(~ DV, space = "free") +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Respondent estimates of Pr(Black | Supporter)",
       x = "",
       y = "")

g1 <- g1 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g1

g2 <- g2 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g2


g <-
  plot_grid(
    g1,
    g2,
    nrow = 2,
    labels = "auto",
    label_fontface = "plain",
    label_fontfamily = "sans",
    label_size = 12
  )
g


g %>%
  ggsave(
    filename = "Fig2.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * (nrow(g1$data))
  )

### ------ Fig. B2 Appendix ------
g1 <-
  gg_df %>%
  filter(DV %in% c("Destruction of property",
                   "Attacks on police officers")) %>%
  ggplot(.,
         aes(
           x = estimate,
           y = Z,
           group = estimator,
           fill = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(position = position_dodgev(height = 0.7),
             size = 3,
             pch = 21) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  scale_color_manual("", values = c("black")) +
  facet_row( ~ DV) +
  scale_x_continuous(labels = percent,
                     expand = c(0, 0),
                     limits = c(0.1, 0.9)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Percentage of respondents that believe most supporters favor this",
       y = "")


g1 <- g1 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g1


g1 %>%
  ggsave(filename = "FigB2.pdf",
         .,
         width = 6.5,
         height = 3)

### ------ Table B2.1 Appendix----------
tab_top <-
  g1$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -term, -.group, -outcome) %>%
  arrange(DV, desc(Z), estimator) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, Z, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry)

tab_top %>%
  write.table(
    .,
    file = "TabB2.1.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

###------- Fig. 3 Manuscript------
the_levs <-
  c(
    "Replacing police with social services",
    "Eliminating police departments",
    "Destruction of property",
    "Attacks on police officers",
    "Reduce use of military equipment",
    "Make firing police easier",
    "Implement Civilian Review Boards",
    "Require police to live where they work",
    "Require Body Worn Cameras",
    "Ban use of chokeholds",
    "Ban use of no-knock warrants",
    "Require warning before shooting"
  )

## Generate plots
g1 <-
  gg_df %>%
  filter(DV %in% the_levs[5:8],
         estimator == "Weighted") %>%
  ggplot(.,
         aes(x = estimate,
             y = Z)) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.7),
    size = 3,
    pch = 21,
    color = "white",
    fill = "black"
  ) +
  facet_col( ~ DV) +
  scale_x_continuous(labels = percent,
                     expand = c(0, 0),
                     limits = c(0.1, 0.99)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "",
       y = "")

g2 <-
  gg_df %>%
  filter(DV %in% the_levs[9:12],
         estimator == "Weighted") %>%
  ggplot(.,
         aes(x = estimate,
             y = Z)) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.7),
    size = 3,
    pch = 21,
    color = "white",
    fill = "black"
  ) +
  facet_col( ~ DV) +
  scale_x_continuous(labels = percent,
                     expand = c(0, 0),
                     limits = c(0.1, 0.99)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.y = element_blank(),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "",
       y = "")

g1 <- g1 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g1

g2 <- g2 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g2

g <-
  plot_grid(
    g1,
    g2,
    ncol = 2,
    label_fontface = "plain",
    label_fontfamily = "sans",
    label_size = 12
  )

g <-
  ggdraw(
    add_sub(
      g,
      "Percentage of respondents that believe most supporters favor this",
      vpadding = grid::unit(0, "lines"),
      y = 6,
      x = 0.5,
      vjust = 4.5,
      size = 10
    )
  )
g

g %>%
  ggsave(filename = "Fig3.pdf",
         .,
         width = 7,
         height = 5)

### ------ Fig. B3 Appendix ------
the_levs <-
  c(
    "Replacing police with social services",
    "Eliminating police departments",
    "Destruction of government property",
    "Attacks on police officers",
    "Reduce use of military equipment",
    "Make firing police easier",
    "Implement Civilian Review Boards",
    "Require police to live where they work",
    "Require Body Worn Cameras",
    "Ban use of chokeholds",
    "Ban use of no-knock warrants",
    "Require warning before shooting"
  )

## Generate plots
g1 <-
  gg_df %>%
  filter(DV %in% the_levs[5:8]) %>%
  ggplot(.,
         aes(
           x = estimate,
           y = Z,
           fill = estimator,
           group = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 2
  ) +
  geom_point(position = position_dodgev(height = 0.7),
             size = 3,
             pch = 21) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  facet_col( ~ DV) +
  scale_x_continuous(
    labels = scales::percent_format(accuracy = 1),
    expand = c(0, 0),
    limits = c(0.1, 0.99)
  ) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "",
       y = "")

g2 <-
  gg_df %>%
  filter(DV %in% the_levs[9:12]) %>%
  ggplot(.,
         aes(
           x = estimate,
           y = Z,
           fill = estimator,
           group = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 2
  ) +
  geom_point(position = position_dodgev(height = 0.7),
             size = 3,
             pch = 21) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  facet_col( ~ DV) +
  scale_x_continuous(
    labels = scales::percent_format(accuracy = 1),
    expand = c(0, 0),
    limits = c(0.1, 0.99)
  ) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.y = element_blank(),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "",
       y = "")

g1 <- g1 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g1

g2 <- g2 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g2

g <-
  plot_grid(
    g1,
    g2,
    ncol = 2,
    labels = NULL,
    label_fontface = "plain",
    label_fontfamily = "sans",
    label_size = 12
  )

g <-
  ggdraw(
    add_sub(
      g,
      "Percentage of respondents that believe most supporters favor this",
      vpadding = grid::unit(0, "lines"),
      y = 6,
      x = 0.5,
      vjust = 4.5,
      size = 10
    )
  )

g

g %>%
  ggsave(filename = "FigB3.pdf",
         .,
         width = 8,
         height = 7)

### ------ Table B3 Appendix ----------
tab_top <-
  g1$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -term, -.group, -outcome) %>%
  arrange(DV, desc(Z), estimator) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, Z, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry)

tab_btm <-
  g2$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -term, -.group, -outcome) %>%
  arrange(DV, desc(Z), estimator) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, Z, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry)

tab_top %>%
  write.table(
    .,
    file = "TabB3.1_top.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

tab_btm %>%
  write.table(
    .,
    file = "TabB3.1_btm.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )


### ------ Fig. 4 Manuscript ----------
g1 <-
  gg_df %>%
  filter(
    DV %in% c(
      "Eliminating police departments",
      "Replacing police with social services"
    ),
    estimator == "Weighted"
  ) %>%
  ggplot(.,
         aes(x = estimate,
             y = Z)) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.7),
    size = 3,
    pch = 21,
    color = "white",
    fill = "black"
  ) +
  facet_row(~ DV) +
  scale_x_continuous(labels = percent,
                     expand = c(0, 0),
                     limits = c(0.1, 0.99)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Percentage of respondents that believe most supporters favor this",
       y = "")


g2 <-
  gg_df %>%
  filter(
    DV %in% c("Spending on police services",
              "Number of police on the street"),
    estimator == "Weighted"
  ) %>%
  ggplot(.,
         aes(x = estimate,
             y = Z)) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.7),
    size = 3,
    pch = 21,
    color = "white",
    fill = "black"
  ) +
  facet_row(~ DV) +
  scale_x_continuous(limits = c(3, 4.5)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x =  "Perceptions of supporters' preferences (1 = Increase a lot,...,5 = Decrease a lot)",
       y = "")

g1 <- g1 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g2 <- g2 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))


g <-
  plot_grid(
    g1,
    g2,
    ncol = 1,
    labels = "auto",
    label_fontface = "plain",
    label_fontfamily = "sans",
    label_size = 12
  )
g

g %>%
  ggsave(
    filename = "Fig4.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * (nrow(g1$data))
  )

### ------ Fig. B4 Appendix----------
g1 <-
  gg_df %>%
  filter(DV %in% c(
    "Eliminating police departments",
    "Replacing police with social services"
  )) %>%
  ggplot(.,
         aes(
           x = estimate,
           y = Z,
           group = estimator,
           fill = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(position = position_dodgev(height = 0.7),
             size = 3,
             pch = 21) +
  facet_row(~ DV) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  scale_x_continuous(labels = percent,
                     expand = c(0, 0),
                     limits = c(0.1, 0.99)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Percentage of respondents that believe most supporters favor this",
       y = "")


g2 <-
  gg_df %>%
  filter(DV %in% c("Spending on police services",
                   "Number of police on the street")) %>%
  ggplot(.,
         aes(
           x = estimate,
           y = Z,
           group = estimator,
           fill = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(position = position_dodgev(height = 0.7),
             size = 3,
             pch = 21) +
  facet_row(~ DV) +
  scale_fill_manual("", values = c("white", "darkgrey"))  +
  scale_x_continuous(limits = c(3, 4.5)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x =  "Perceptions of supporters' preferences (1 = Increase a lot,...,5 = Decrease a lot)",
       y = "")

g1 <- g1 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g2 <- g2 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))


g <-
  plot_grid(
    g1,
    g2,
    ncol = 1,
    labels = "auto",
    label_fontface = "plain",
    label_fontfamily = "sans",
    label_size = 12
  )
g

g %>%
  ggsave(
    filename = "FigB4.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * (nrow(g1$data))
  )

### ------ Table B4 Appendix------
tab_top <-
  g1$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -term, -.group, -outcome) %>%
  arrange(DV, desc(Z), estimator) %>%
  mutate(entry = table_entry(est = estimate, se = std.error, digits = 2)) %>%
  select(DV, estimator, Z, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry)

tab_btm <-
  g2$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -term, -.group, -outcome) %>%
  arrange(DV, desc(Z), estimator) %>%
  mutate(entry = table_entry(est = estimate, se = std.error, digits = 2)) %>%
  select(DV, estimator, Z, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry)

tab_top %>%
  write.table(
    .,
    file = "TabB4.1_top.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

tab_btm %>%
  write.table(
    .,
    file = "TabB4.1_btm.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

### ------ Fig. B6 Appendix ------
these <-
  c(
    "Responding to violent crimes",
    "Responding to property crimes",
    "Responding to vice crimes",
    "Crime prevention",
    "Maintaining civil order",
    "Conducting vehicle stops",
    "Providing social services",
    "Community engagement"
  )

g1 <-
  gg_df %>%
  filter(DV %in% these) %>%
  mutate(DV = factor(DV, levels = these)) %>%
  ggplot(.,
         aes(
           x = estimate,
           y = Z,
           fill = estimator,
           group = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.8),
    size = 2
  ) +
  geom_point(position = position_dodgev(height = 0.8),
             size = 3,
             pch = 21) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  facet_wrap(~ DV, ncol = 2) +
  scale_x_continuous(labels = percent,
                     expand = c(0, 0),
                     limits = c(0.1, 0.99)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10),
    plot.margin = grid::unit(c(0, 0, 0, 0), "mm")
  ) +
  labs(subtitle = "",
       x = "Percentage that believe most supporters oppose police doing this",
       y = "")


g1 %>%
  ggsave(filename = "FigB6.pdf",
         .,
         width = 7,
         height = 6)

### ------ Table B6 Appendix----------
g1$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -term, -.group, -outcome) %>%
  arrange(DV, desc(Z), estimator) %>%
  head(., 24) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, Z, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry) %>%
  write.table(
    .,
    file = "TabB6.1_top.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

g1$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -term, -.group, -outcome) %>%
  arrange(DV, desc(Z), estimator) %>%
  tail(., 24) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, Z, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry) %>%
  write.table(
    .,
    file = "TabB6.1_btm.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )


### ------ Fig. B7 Appendix -------

## Now split out categories for bottom
call_df <-
  dat %>%
  select(starts_with("call_matrix"), weights, dem_pid_3) %>%
  pivot_longer(
    cols = starts_with("call_matrix"),
    values_to = "Y",
    names_to = "DV"
  ) %>%
  mutate(DV = str_replace(DV, "call_matrix_", ""),
         DV = factor(DV, levels = rev(
           c(
             "domestic",
             "fight",
             "suspicious",
             "counterfit",
             "racist",
             "drunk",
             "heroin",
             "teens"
           )
         ),
         labels = rev(
           c(
             "Person threatening spouse with knife",
             "Fist fight outside bar",
             "Person looking in car windows",
             "Person using counterfit money",
             "Person yelling racist slurs at voting queue",
             "Drunk person yelling at outdoor diners",
             "Teens playing loud music at night in park",
             "Heroin overdose in park"
           )
         )))

call_nothing <-
  call_df %>%
  filter(!is.na(Y), !is.na(DV)) %>%
  group_by(DV) %>%
  do(tidy(lm_robust(Y == "Do nothing" ~ 1, data = .))) %>%
  mutate(outcome = "Do nothing")


call_police_plus <-
  call_df %>%
  filter(!is.na(Y), !is.na(DV)) %>%
  group_by(DV) %>%
  do(tidy(lm_robust(
    Y == "Send the police and someone else" ~ 1, data = .
  ))) %>%
  mutate(outcome = "Send both")


call_nonpolice <-
  call_df %>%
  filter(!is.na(Y), !is.na(DV)) %>%
  group_by(DV) %>%
  do(tidy(
    lm_robust(
      Y == "Send someone other than the police, like the fire department or social services" ~ 1,
      data = .
    )
  )) %>%
  mutate(outcome = "Send someone else")

call_police <-
  call_df %>%
  filter(!is.na(Y), !is.na(DV)) %>%
  group_by(DV) %>%
  do(tidy(lm_robust(Y == "Send the police" ~ 1, data = .))) %>%
  mutate(outcome = "Send the police")


gg_u_df <-
  bind_rows(call_nothing, call_police_plus, call_police, call_nonpolice) %>%
  mutate(estimator = "Unweighted") %>%
  select(-term)

## Repeat w/ weights
call_nothing <-
  call_df %>%
  filter(!is.na(Y), !is.na(DV)) %>%
  group_by(DV) %>%
  do(tidy(lm_robust(
    Y == "Do nothing" ~ 1, data = .,
    weights = weights
  ))) %>%
  mutate(outcome = "Do nothing")


call_police_plus <-
  call_df %>%
  filter(!is.na(Y), !is.na(DV)) %>%
  group_by(DV) %>%
  do(tidy(
    lm_robust(
      Y == "Send the police and someone else" ~ 1,
      data = .,
      weights = weights
    )
  )) %>%
  mutate(outcome = "Send both")


call_nonpolice <-
  call_df %>%
  filter(!is.na(Y), !is.na(DV)) %>%
  group_by(DV) %>%
  do(tidy(
    lm_robust(
      Y == "Send someone other than the police, like the fire department or social services" ~ 1,
      data = .
    )
  )) %>%
  mutate(outcome = "Send someone else")

call_police <-
  call_df %>%
  filter(!is.na(Y), !is.na(DV)) %>%
  group_by(DV) %>%
  do(tidy(lm_robust(
    Y == "Send the police" ~ 1, data = ., weights = weights
  ))) %>%
  mutate(outcome = "Send the police")


gg_w_df <-
  bind_rows(call_nothing, call_police_plus, call_police, call_nonpolice) %>%
  mutate(estimator = "Weighted") %>%
  select(-term)


## Combine and plot
the_levs <-
  c(
    "Person threatening spouse with knife",
    "Fist fight outside bar",
    "Person looking in car windows",
    "Person using counterfit money",
    "Person yelling racist slurs at voting queue",
    "Drunk person yelling at outdoor diners",
    "Teens playing loud music at night in park",
    "Heroin overdose in park"
  )

the_labs <-
  c(
    "Person threatening spouse with knife",
    "Fist fight outside bar",
    "Person looking in car windows",
    "Person using counterfit money",
    "Person yelling racist slurs at voters",
    "Drunk person yelling at diners",
    "Teens playing loud music in park",
    "Heroin overdose in park"
  )

gg_btm_split <-
  bind_rows(gg_u_df, gg_w_df) %>%
  mutate(activity = factor(DV, levels = the_levs,
                           labels = the_labs))

## Sanity check
gg_btm_split %>%
  group_by(DV, estimator) %>%
  summarise(check = sum(estimate))


g3 <-
  gg_btm_split %>%
  ggplot(.,
         aes(
           x = estimate,
           y = reorder(outcome, estimate, na.rm = TRUE),
           fill = estimator,
           group = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 1),
    size = 2
  ) +
  geom_point(position = position_dodgev(height = 1),
             size = 3,
             pch = 21) +
  facet_wrap(~ activity, ncol = 3, scales = "free_x") +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  scale_x_continuous(
    labels = percent_format(accuracy = 1),
    expand = c(0, 0),
    limits = c(0, 0.9)
  ) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Percentage of respondents supporting each response type",
       y = "")

g3 %>%
  ggsave(filename = "FigB7.pdf",
         .,
         width = 9,
         height = 6)

### ------ Table B7.1 Appendix ----------
g3$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -.group, -activity) %>%
  arrange(desc(DV), estimator, desc(outcome)) %>%
  head(., 32) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, outcome, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry) %>%
  write.table(
    .,
    file = "TabB7_top.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

g3$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -.group, -activity) %>%
  arrange(desc(DV), estimator, desc(outcome)) %>%
  tail(., 32) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, outcome, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry) %>%
  write.table(
    .,
    file = "TabB7_btm.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

### ------ Fig. B9 Appendix-----
these <- c(
  "Spending on education",
  "Spending on healthcare",
  "Spending on housing",
  "Spending on social services"
)

g <-
  gg_df %>%
  filter(DV %in% these) %>%
  ggplot(.,
         aes(
           x = estimate,
           y = Z,
           fill = estimator,
           group = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.9),
    size = 2
  ) +
  geom_point(position = position_dodgev(height = 0.9),
             size = 3,
             pch = 21) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  facet_col( ~ DV) +
  scale_x_continuous(limits = c(1.5, 3)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_blank()
  ) +
  labs(subtitle = "",
       x = "Perceptions of supporters' preferences\n(1 = Increase a lot,...,5 = Decrease a lot)",
       y = "")

g <- g + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g


g %>%
  ggsave(filename = "FigB9.pdf",
         .,
         width = 6.5,
         height = 7)

### ------ Table B9.1 Appendix----------
g$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -.group, -outcome) %>%
  arrange(DV, desc(Z), estimator) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, Z, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry) %>%
  write.table(
    .,
    file = "TabB9.1.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

### ------ Estimate treatment effects ----------
## Difference in means.
mod_dim_u <-
  dat %>%
  select(
    Z,
    othr_supp_eliminate_n,
    othr_supp_replace_n,
    othr_spend_police_n,
    othr_change_numcops_n,
    othr_supp_destroy_n,
    othr_supp_attack_n,
    othr_supp_verbal_n,
    othr_supp_choke_n,
    othr_supp_knock_n,
    othr_supp_military_n,
    othr_supp_identify_n,
    othr_supp_crb_n,
    othr_supp_bwc_n,
    othr_supp_live_n,
    othr_supp_protest_brutal_n,
    othr_supp_protest_blm_n,
    othr_supp_protest_defund_n,
    othr_act_violent_n,
    othr_act_property_n,
    othr_act_vice_n,
    othr_act_order_n,
    othr_act_prevent_n,
    othr_act_traffic_n,
    othr_act_social_n,
    othr_act_engage_n,
    othr_spend_edu_n,
    othr_spend_health_n,
    othr_spend_housing_n,
    othr_spend_social_n
  ) %>%
  pivot_longer(
    othr_supp_eliminate_n:othr_spend_social_n,
    names_to = "DV",
    values_to = "Y"
  ) %>%
  group_by(DV) %>%
  do(tidy(lm_robust(Y ~ Z, data = .)))



mod_dim_w <-
  dat %>%
  select(
    Z,
    weights,
    othr_supp_eliminate_n,
    othr_supp_replace_n,
    othr_spend_police_n,
    othr_change_numcops_n,
    othr_supp_destroy_n,
    othr_supp_attack_n,
    othr_supp_verbal_n,
    othr_supp_choke_n,
    othr_supp_knock_n,
    othr_supp_military_n,
    othr_supp_identify_n,
    othr_supp_crb_n,
    othr_supp_bwc_n,
    othr_supp_live_n,
    othr_supp_protest_brutal_n,
    othr_supp_protest_blm_n,
    othr_supp_protest_defund_n,
    othr_act_violent_n,
    othr_act_property_n,
    othr_act_vice_n,
    othr_act_order_n,
    othr_act_prevent_n,
    othr_act_traffic_n,
    othr_act_social_n,
    othr_act_engage_n,
    othr_spend_edu_n,
    othr_spend_health_n,
    othr_spend_housing_n,
    othr_spend_social_n
  ) %>%
  pivot_longer(
    othr_supp_eliminate_n:othr_spend_social_n,
    names_to = "DV",
    values_to = "Y"
  ) %>%
  group_by(DV) %>%
  do(tidy(lm_robust(Y ~ Z, data = ., weights = weights)))


## Covariate-adjusted
mod_adj_u <-
  dat %>%
  select(
    covariates,
    Z,
    othr_supp_eliminate_n,
    othr_supp_replace_n,
    othr_spend_police_n,
    othr_change_numcops_n,
    othr_supp_destroy_n,
    othr_supp_attack_n,
    othr_supp_verbal_n,
    othr_supp_choke_n,
    othr_supp_knock_n,
    othr_supp_military_n,
    othr_supp_identify_n,
    othr_supp_crb_n,
    othr_supp_bwc_n,
    othr_supp_live_n,
    othr_supp_protest_brutal_n,
    othr_supp_protest_blm_n,
    othr_supp_protest_defund_n,
    othr_act_violent_n,
    othr_act_property_n,
    othr_act_vice_n,
    othr_act_order_n,
    othr_act_prevent_n,
    othr_act_traffic_n,
    othr_act_social_n,
    othr_act_engage_n,
    othr_spend_edu_n,
    othr_spend_health_n,
    othr_spend_housing_n,
    othr_spend_social_n
  ) %>%
  pivot_longer(
    othr_supp_eliminate_n:othr_spend_social_n,
    names_to = "DV",
    values_to = "Y"
  ) %>%
  group_by(DV) %>%
  do(tidy(lm_lin(
    Y ~ Z, data = .,
    covariates = as.formula(paste("~", paste(
      covariates, collapse = "+"
    )))
  )))

mod_adj_w <-
  dat %>%
  select(
    covariates,
    Z,
    weights,
    othr_supp_eliminate_n,
    othr_supp_replace_n,
    othr_spend_police_n,
    othr_change_numcops_n,
    othr_supp_destroy_n,
    othr_supp_attack_n,
    othr_supp_verbal_n,
    othr_supp_choke_n,
    othr_supp_knock_n,
    othr_supp_military_n,
    othr_supp_identify_n,
    othr_supp_crb_n,
    othr_supp_bwc_n,
    othr_supp_live_n,
    othr_supp_protest_brutal_n,
    othr_supp_protest_blm_n,
    othr_supp_protest_defund_n,
    othr_act_violent_n,
    othr_act_property_n,
    othr_act_vice_n,
    othr_act_order_n,
    othr_act_prevent_n,
    othr_act_traffic_n,
    othr_act_social_n,
    othr_act_engage_n,
    othr_spend_edu_n,
    othr_spend_health_n,
    othr_spend_housing_n,
    othr_spend_social_n
  ) %>%
  pivot_longer(
    othr_supp_eliminate_n:othr_spend_social_n,
    names_to = "DV",
    values_to = "Y"
  ) %>%
  group_by(DV) %>%
  do(tidy(lm_lin(
    Y ~ Z,
    data = .,
    weights = weights,
    covariates = as.formula(paste("~", paste(
      covariates, collapse = "+"
    )))
  )))

the_levs <- c(
  "othr_supp_destroy_n",
  "othr_supp_attack_n",
  "othr_supp_eliminate_n",
  "othr_supp_replace_n",
  "othr_spend_police_n",
  "othr_change_numcops_n",
  "othr_supp_military_n",
  "othr_supp_identify_n",
  "othr_supp_crb_n",
  "othr_supp_live_n",
  "othr_supp_bwc_n",
  "othr_supp_choke_n",
  "othr_supp_knock_n",
  "othr_supp_verbal_n",
  "othr_supp_protest_brutal_n",
  "othr_supp_protest_blm_n",
  "othr_supp_protest_defund_n",
  "othr_act_violent_n",
  "othr_act_property_n",
  "othr_act_vice_n",
  "othr_act_order_n",
  "othr_act_prevent_n",
  "othr_act_traffic_n",
  "othr_act_social_n",
  "othr_act_engage_n",
  "othr_spend_edu_n",
  "othr_spend_health_n",
  "othr_spend_housing_n",
  "othr_spend_social_n"
)

the_labs =  c(
  "Destruction of property",
  "Attacks on police officers",
  "Eliminating police departments",
  "Replacing police with social services",
  "Spending on police services",
  "Number of police on the street",
  "Reduce use of military equipment",
  "Make firing police easier",
  "Implement Civilian Review Boards",
  "Require police to live where they work",
  "Require Body Worn Cameras",
  "Ban use of chokeholds",
  "Ban use of no-knock warrants",
  "Require warning before shooting",
  "Protest against Police Brutality",
  "Protests by Black Lives Matter",
  "Protests to Defund the Police",
  "Responding to violent crimes",
  "Responding to property crimes",
  "Responding to vice crimes",
  "Conducting vehicle stops",
  "Crime prevention",
  "Maintaining civil order",
  "Providing social services",
  "Community engagement",
  "Spending on education",
  "Spending on healthcare",
  "Spending on housing",
  "Spending on social services"
)

est_df <-
  bind_rows(
    mod_dim_u %>%
      mutate(estimator = "Diff-in-Means",
             weighted = FALSE),
    
    mod_dim_w %>%
      mutate(estimator = "Diff-in-Means",
             weighted = TRUE),
    
    mod_adj_u %>%
      mutate(estimator = "Covariate-adjusted",
             weighted = FALSE),
    
    mod_adj_w %>%
      mutate(estimator = "Covariate-adjusted",
             weighted = TRUE)
  ) %>%
  select(-df) %>%
  filter(term %in% c("ZDefund", "ZAbolish", "(Intercept)")) %>%
  mutate(
    term = case_when(
      term == "(Intercept)" ~ "Reform Mean",
      term == "ZDefund" ~ "Defund",
      term == "ZAbolish" ~ "Abolish"
    ),
    outcome = factor(DV, levels = the_levs, labels = the_labs)
  )

## Generate tables for word
ols_table <-
  est_df %>%
  filter(term != "Reform Mean")

### ------ Table B2.2 Appendix ----------
ols_table %>%
  filter(outcome %in% the_labs[1:2]) %>%
  arrange(desc(estimator), outcome, term, weighted) %>%
  group_by(estimator, weighted) %>%
  mutate(entry = make_entry(
    est = estimate,
    se = std.error,
    p = p.adjust(p.value, "fdr"),
    digits = 2
  )) %>%
  dplyr::select(term, estimator, outcome, weighted, entry) %>%
  pivot_wider(names_from = outcome:weighted,
              values_from = c(entry)) %>%
  write.table(
    .,
    file = "TabB2.2.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )


### ------ Table B3.2 Appendix ----------
tmp <-
  ols_table %>%
  filter(outcome %in% the_labs[7:14]) %>%
  arrange(desc(estimator), outcome, term, weighted) %>%
  group_by(estimator, weighted) %>%
  mutate(entry = make_entry(
    est = estimate,
    se = std.error,
    p = p.adjust(p.value, "fdr"),
    digits = 2
  ))

tmp %>%
  filter(outcome %in% the_labs[7:10]) %>%
  dplyr::select(term, estimator, outcome, weighted, entry) %>%
  pivot_wider(names_from = outcome:weighted,
              values_from = c(entry)) %>%
  write.table(
    .,
    file = "TabB3.2_top.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

tmp %>%
  filter(outcome %in% the_labs[11:14]) %>%
  dplyr::select(term, estimator, outcome, weighted, entry) %>%
  pivot_wider(names_from = outcome:weighted,
              values_from = c(entry)) %>%
  write.table(
    .,
    file = "TabB3.2_btm.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

### ------ Table B4.2 Appendix ----------
ols_table %>%
  filter(outcome %in% the_labs[3:4]) %>%
  arrange(desc(estimator), outcome, term, weighted) %>%
  group_by(estimator, weighted) %>%
  mutate(entry = make_entry(
    est = estimate,
    se = std.error,
    p = p.adjust(p.value, "fdr"),
    digits = 2
  )) %>%
  dplyr::select(term, estimator, outcome, weighted, entry) %>%
  pivot_wider(names_from = outcome:weighted,
              values_from = c(entry)) %>%
  write.table(
    .,
    file = "TabB4.2_top.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )


ols_table %>%
  filter(outcome %in% the_labs[5:6]) %>%
  arrange(desc(estimator), outcome, term, weighted) %>%
  group_by(estimator, weighted) %>%
  mutate(entry = make_entry(
    est = estimate,
    se = std.error,
    p = p.adjust(p.value, "fdr"),
    digits = 2
  )) %>%
  dplyr::select(term, estimator, outcome, weighted, entry) %>%
  pivot_wider(names_from = outcome:weighted,
              values_from = c(entry)) %>%
  write.table(
    .,
    file = "TabB4.2_btm.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

### ------ Table B6.2 Appendix ----------
ols_table %>%
  filter(outcome %in% the_labs[c(18:20, 22)]) %>%
  arrange(desc(estimator), outcome, term, weighted) %>%
  group_by(estimator, weighted) %>%
  mutate(entry = make_entry(
    est = estimate,
    se = std.error,
    p = p.adjust(p.value, "fdr"),
    digits = 2
  )) %>%
  dplyr::select(term, estimator, outcome, weighted, entry) %>%
  pivot_wider(names_from = outcome:weighted,
              values_from = c(entry)) %>%
  write.table(
    .,
    file = "TabB6.2_top.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

ols_table %>%
  filter(outcome %in% the_labs[c(23, 21, 24:25)]) %>%
  arrange(desc(estimator), outcome, term, weighted) %>%
  group_by(estimator, weighted) %>%
  mutate(entry = make_entry(
    est = estimate,
    se = std.error,
    p = p.adjust(p.value, "fdr"),
    digits = 2
  )) %>%
  dplyr::select(term, estimator, outcome, weighted, entry) %>%
  pivot_wider(names_from = outcome:weighted,
              values_from = c(entry)) %>%
  write.table(
    .,
    file = "TabB6.2_btm.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

### ------ Table B9.2 Appendix ----------
ols_table %>%
  filter(outcome %in% the_labs[26:29]) %>%
  arrange(desc(estimator), outcome, term, weighted) %>%
  group_by(estimator, weighted) %>%
  mutate(entry = make_entry(
    est = estimate,
    se = std.error,
    p = p.adjust(p.value, "fdr"),
    digits = 2
  )) %>%
  dplyr::select(term, estimator, outcome, weighted, entry) %>%
  pivot_wider(names_from = outcome:weighted,
              values_from = c(entry)) %>%
  write.table(
    .,
    file = "TabB9.2.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

### ------ Fig. 5 Manuscript-----
call_df <-
  dat %>%
  select(starts_with("call_matrix"), weights, dem_pid_3) %>%
  pivot_longer(
    cols = starts_with("call_matrix"),
    values_to = "Y",
    names_to = "DV"
  ) %>%
  mutate(
    DV = str_replace(DV, "call_matrix_", ""),
    DV = factor(DV, levels = rev(
      c(
        "domestic",
        "fight",
        "suspicious",
        "counterfit",
        "racist",
        "drunk",
        "heroin",
        "teens"
      )
    ),
    labels = rev(
      c(
        "Person threatening spouse with knife",
        "Fist fight outside bar",
        "Person looking in car windows",
        "Person using counterfit money",
        "Person yelling racist slurs at voting queue",
        "Drunk person yelling at outdoor diners",
        "Teens playing loud music at night in park",
        "Heroin overdose in park"
      )
    )),
    Y = case_when(
      Y == "Do nothing" ~ "No police involvement",
      Y ==  "Send someone other than the police, like the fire department or social services" ~ "No police involvement",
      Y == "Send the police and someone else" ~ "Police involvement",
      Y == "Send the police" ~ "Police involvement"
    ),
    Y = factor(Y,
               levels = rev(
                 c("No police involvement",
                   "Police involvement")
               ))
  )

call_police <-
  call_df %>%
  filter(!is.na(Y), !is.na(DV)) %>%
  group_by(DV) %>%
  do(tidy(lm_robust(Y == "Police involvement" ~ 1, data = .))) %>%
  mutate(outcome = "Police involvement")

call_nonpolice <-
  call_df %>%
  filter(!is.na(Y), !is.na(DV)) %>%
  group_by(DV) %>%
  do(tidy(lm_robust(Y == "No police involvement" ~ 1, data = .))) %>%
  mutate(outcome = "No police involvement")


gg_u_df <-
  bind_rows(call_police, call_nonpolice) %>%
  mutate(estimator = "Unweighted") %>%
  select(-term)

## Repeat w/ weights
call_police <-
  call_df %>%
  filter(!is.na(Y), !is.na(DV)) %>%
  group_by(DV) %>%
  do(tidy(lm_robust(
    Y == "Police involvement" ~ 1,
    data = .,
    weights = weights
  ))) %>%
  mutate(outcome = "Police involvement")

call_nonpolice <-
  call_df %>%
  filter(!is.na(Y), !is.na(DV)) %>%
  group_by(DV) %>%
  do(tidy(
    lm_robust(
      Y == "No police involvement" ~ 1,
      data = .,
      weights = weights
    )
  )) %>%
  mutate(outcome = "No police involvement")


gg_w_df <-
  bind_rows(call_police, call_nonpolice) %>%
  mutate(estimator = "Weighted") %>%
  select(-term)

## Combine and plot
the_levs <-
  c(
    "Person threatening spouse with knife",
    "Fist fight outside bar",
    "Person looking in car windows",
    "Person using counterfit money",
    "Person yelling racist slurs at voting queue",
    "Drunk person yelling at outdoor diners",
    "Teens playing loud music at night in park",
    "Heroin overdose in park"
  )

the_labs <-
  c(
    "Person threatening spouse with knife",
    "Fist fight outside bar",
    "Person looking in car windows",
    "Person using counterfit money",
    "Person yelling racist slurs at voting queue",
    "Drunk person yelling at outdoor diners",
    "Teens playing loud music at night in park",
    "Heroin overdose in park"
  )

gg_top <-
  bind_rows(gg_u_df, gg_w_df) %>%
  mutate(activity = factor(DV, levels = rev(the_levs),
                           labels = rev(the_labs)))

## Sanity check
gg_top %>%
  group_by(DV, estimator) %>%
  summarise(check = sum(estimate))


## Plot results:
g1 <-
  gg_top %>%
  filter(outcome == "No police involvement", estimator == "Weighted") %>%
  ggplot(.,
         aes(
           x = estimate,
           y = activity,
           fill = estimator,
           group = estimator,
           color = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.7),
    size = 3,
    pch = 21,
    color = "white"
  ) +
  scale_fill_manual("", values = c("black")) +
  scale_color_manual("", values = c("black")) +
  scale_x_continuous(
    labels = percent_format(accuracy = 1),
    expand = c(0, 0),
    limits = c(0, 0.55)
  ) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Percentage that believe police should not be involved",
       y = "")

## Bottom of plot
gg_u_df <-
  dat %>%
  select(Z_cut,
         weights,
         cuts_crime_binary,
         cuts_safe_binary) %>%
  pivot_longer(cuts_crime_binary:cuts_safe_binary,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(Z_cut, DV) %>%
  do(tidy(lm_robust(Y ~ 1, data = .))) %>%
  filter(term == "(Intercept)") %>%
  ungroup() %>%
  mutate(estimator = "Unweighted")

gg_w_df <-
  dat %>%
  select(Z_cut,
         weights,
         cuts_crime_binary,
         cuts_safe_binary) %>%
  pivot_longer(cuts_crime_binary:cuts_safe_binary,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(Z_cut, DV) %>%
  do(tidy(lm_robust(Y ~ 1, data = ., weights = weights))) %>%
  filter(term == "(Intercept)") %>%
  ungroup() %>%
  mutate(estimator = "Weighted")

gg_btm <-
  bind_rows(gg_u_df, gg_w_df) %>%
  mutate(
    Z_cut = factor(Z_cut,
                   levels = rev(c(1, 2, 3, 4)),
                   labels = rev(
                     c(
                       "10% reduction",
                       "25% reduction",
                       "50% reduction",
                       "100% reduction"
                     )
                   )),
    DV = factor(
      DV,
      levels = c("cuts_crime_binary", "cuts_safe_binary"),
      labels = c("Increased crime", "Decreased safety")
    )
  )

g2 <-
  gg_btm %>%
  filter(estimator == "Weighted") %>%
  ggplot(.,
         aes(
           x = estimate,
           y = Z_cut,
           fill = estimator,
           group = estimator,
           color = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.7),
    size = 3,
    color = "white",
    pch = 21
  ) +
  scale_fill_manual("", values = c("black")) +
  scale_color_manual("", values = c("black")) +
  facet_row( ~ DV) +
  scale_x_continuous(labels = percent_format(accuracy = 1),
                     limits = c(0.45, 0.86)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Percentage that believe reducing the number of police will cause this",
       y = "")


g1 <- g1 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g2 <- g2 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))


g <-
  plot_grid(
    g1,
    g2,
    nrow = 2,
    labels = "auto",
    label_fontface = "plain",
    label_fontfamily = "sans",
    label_size = 12,
    rel_heights = c(1.25, 1)
  )
g


g %>%
  ggsave(
    filename = "Fig5.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * nrow(g1$data)
  )

### ------ Fig. B5 Appendix----------
g1 <-
  gg_top %>%
  filter(outcome == "No police involvement") %>%
  ggplot(.,
         aes(
           x = estimate,
           y = activity,
           fill = estimator,
           group = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(position = position_dodgev(height = 0.7),
             size = 3,
             pch = 21) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  scale_x_continuous(
    labels = percent_format(accuracy = 1),
    expand = c(0, 0),
    limits = c(0, 0.55)
  ) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Percentage that believe police should not be involved",
       y = "")

g2 <-
  gg_btm %>%
  ggplot(.,
         aes(
           x = estimate,
           y = Z_cut,
           fill = estimator,
           group = estimator
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(position = position_dodgev(height = 0.7),
             size = 3,
             pch = 21) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  facet_row( ~ DV) +
  scale_x_continuous(labels = percent_format(accuracy = 1),
                     limits = c(0.45, 0.86)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Percentage that believe reducing the number of police will cause this",
       y = "")


g1 <- g1 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g2 <- g2 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))


g <-
  plot_grid(
    g1,
    g2,
    nrow = 2,
    labels = "auto",
    label_fontface = "plain",
    label_fontfamily = "sans",
    label_size = 12,
    rel_heights = c(1.25, 1)
  )
g


g %>%
  ggsave(
    filename = "FigB5.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * nrow(g1$data)
  )

### ------ Table B5 Appendix ------
g1$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -.group, -activity, -outcome) %>%
  arrange(desc(DV), estimator) %>%
  head(., 8) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry) %>%
  write.table(
    .,
    file = "TabB5.1a_top.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

g1$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -.group, -activity, -outcome) %>%
  arrange(desc(DV), estimator) %>%
  tail(., 8) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry) %>%
  write.table(
    .,
    file = "TabB5.1a_btm.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

g2$data %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -outcome) %>%
  arrange(DV, desc(Z_cut), estimator) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, Z_cut, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry) %>%
  write.table(
    .,
    file = "TabB5.1b_top.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

## Estimates for DV without binary conversion
gg_u_df <-
  dat %>%
  select(Z_cut,
         weights,
         dem_pid_7n,
         cuts_crime_n,
         cuts_safe_n,
         dem_race_4) %>%
  pivot_longer(cuts_crime_n:cuts_safe_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(Z_cut, DV) %>%
  do(tidy(lm_robust(Y ~ 1, data = .))) %>%
  filter(term == "(Intercept)") %>%
  ungroup() %>%
  mutate(estimator = "Unweighted")

gg_w_df <-
  dat %>%
  select(Z_cut,
         weights,
         dem_pid_7n,
         cuts_crime_n,
         cuts_safe_n,
         dem_race_4) %>%
  pivot_longer(cuts_crime_n:cuts_safe_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(Z_cut, DV) %>%
  do(tidy(lm_robust(Y ~ 1, data = ., weights = weights))) %>%
  filter(term == "(Intercept)") %>%
  ungroup() %>%
  mutate(estimator = "Weighted")

gg_btm <-
  bind_rows(gg_u_df, gg_w_df) %>%
  mutate(
    Z_cut = factor(Z_cut,
                   levels = rev(c(1, 2, 3, 4)),
                   labels = rev(
                     c(
                       "10% reduction",
                       "25% reduction",
                       "50% reduction",
                       "100% reduction"
                     )
                   )),
    DV = factor(
      DV,
      levels = c("cuts_crime_n", "cuts_safe_n"),
      labels = c("Increased crime", "Decreased safety")
    )
  )


gg_btm %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -outcome) %>%
  arrange(DV, desc(Z_cut), estimator) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, estimator, Z_cut, entry) %>%
  pivot_wider(names_from = DV:estimator,
              values_from = entry) %>%
  write.table(
    .,
    file = "TabB5.1b_btm.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

### ------ Table B5.2 Appendix ----------

## Difference in means.
mod_dim_u <-
  dat %>%
  select(Z_cut,
         cuts_safe_binary,
         cuts_crime_binary,
         cuts_safe_n,
         cuts_crime_n,) %>%
  pivot_longer(cuts_safe_binary:cuts_crime_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(DV) %>%
  do(tidy(lm_robust(Y ~ Z_cut, data = .)))

mod_dim_w <-
  dat %>%
  select(Z_cut,
         weights,
         cuts_safe_binary,
         cuts_crime_binary,
         cuts_safe_n,
         cuts_crime_n,) %>%
  pivot_longer(cuts_safe_binary:cuts_crime_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(DV) %>%
  do(tidy(lm_robust(
    Y ~ Z_cut, weights = weights, data = .
  )))

## Covariate-adjusted
mod_adj_u <-
  dat %>%
  select(covariates,
         Z_cut,
         cuts_safe_binary,
         cuts_crime_binary,
         cuts_safe_n,
         cuts_crime_n,) %>%
  pivot_longer(cuts_safe_binary:cuts_crime_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(DV) %>%
  do(tidy(lm_lin(
    Y ~ Z_cut, data = .,
    covariates = as.formula(paste("~", paste(
      covariates, collapse = "+"
    )))
  )))

mod_adj_w <-
  dat %>%
  select(
    covariates,
    Z_cut,
    weights,
    cuts_safe_binary,
    cuts_crime_binary,
    cuts_safe_n,
    cuts_crime_n,
  ) %>%
  pivot_longer(cuts_safe_binary:cuts_crime_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(DV) %>%
  do(tidy(
    lm_lin(
      Y ~ Z_cut,
      data = .,
      weights = weights,
      covariates = as.formula(paste("~", paste(
        covariates, collapse = "+"
      )))
    )
  ))

the_levs <- c("cuts_safe_binary",
              "cuts_crime_binary",
              "cuts_safe_n",
              "cuts_crime_n")

the_labs =  c("Increased crime",
              "Decreased safety",
              "Crime levels",
              "Safety levels")

est_cuts_df <-
  bind_rows(
    mod_dim_u %>%
      mutate(estimator = "Diff-in-Means",
             weighted = FALSE),
    
    mod_dim_w %>%
      mutate(estimator = "Diff-in-Means",
             weighted = TRUE),
    
    mod_adj_u %>%
      mutate(estimator = "Covariate-adjusted",
             weighted = FALSE),
    
    mod_adj_w %>%
      mutate(estimator = "Covariate-adjusted",
             weighted = TRUE)
  ) %>%
  select(-df) %>%
  filter(term %in% c("Z_cut25", "Z_cut50", "Z_cut100", "(Intercept)")) %>%
  mutate(
    term = case_when(
      term == "(Intercept)" ~ "10% Mean",
      term == "Z_cut25" ~ "25%",
      term == "Z_cut50" ~ "50%",
      term == "Z_cut100" ~ "100%"
    ),
    outcome = factor(DV, levels = the_levs, labels = the_labs)
  )


## Generate tables for word
ols_table <-
  est_cuts_df %>%
  filter(term != "10% Mean")

## Table B5.2 top
ols_table %>%
  filter(outcome %in% c(the_labs[1:2])) %>%
  arrange(desc(estimator), outcome, term, weighted) %>%
  group_by(estimator, weighted) %>%
  mutate(entry = make_entry(
    est = estimate,
    se = std.error,
    p = p.adjust(p.value, "fdr"),
    digits = 2
  )) %>%
  dplyr::select(term, estimator, outcome, weighted, entry) %>%
  pivot_wider(names_from = outcome:weighted,
              values_from = c(entry)) %>%
  write.table(
    .,
    file = "TabB5.2_top.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

## Table B5.2 bottom
ols_table %>%
  filter(outcome %in% c(the_labs[2:4])) %>%
  arrange(desc(estimator), outcome, term, weighted) %>%
  group_by(estimator, weighted) %>%
  mutate(entry = make_entry(
    est = estimate,
    se = std.error,
    p = p.adjust(p.value, "fdr"),
    digits = 2
  )) %>%
  dplyr::select(term, estimator, outcome, weighted, entry) %>%
  pivot_wider(names_from = outcome:weighted,
              values_from = c(entry)) %>%
  write.table(
    .,
    file = "TabB5.2_btm.txt",
    sep = ",",
    quote = FALSE,
    row.names = F,
    col.names = FALSE
  )

### ------ Fig. B8.1 Appendix ------
gg_df <-
  dat %>%
  select(Z,
         weights,
         othr_spend_police_n) %>%
  rename(DV = othr_spend_police_n) %>%
  mutate(
    DV = case_when(
      DV == 5 ~ "Increase a lot",
      DV == 4 ~ "Increase a little",
      DV == 3 ~ "Keep the same",
      DV == 2 ~ "Decrease a little",
      DV == 1 ~ "Decrease a lot"
    ),
    DV = factor(
      DV,
      levels = c(
        "Increase a lot",
        "Increase a little",
        "Keep the same",
        "Decrease a little",
        "Decrease a lot"
      )
    ),
    Z = factor(
      Z,
      levels = c("Abolish", "Defund", "Reform"),
      labels = c("Abolishing", "Defunding", "Reforming")
    )
  ) %>%
  filter(!is.na(DV)) %>%
  group_by(Z, DV) %>%
  summarise(n_wt = sum(weights),
            n = n()) %>%
  group_by(Z) %>%
  mutate(weighted = n_wt / sum(n_wt),
         unweighted = n / sum(n)) %>%
  pivot_longer(cols = weighted:unweighted,
               names_to = "type",
               values_to = "prop")

## Sanity check
gg_df %>% group_by(Z, type) %>% summarise(sum(prop))

g1 <-
  gg_df %>%
  ggplot(., aes(y = prop, x = DV, fill = type)) +
  geom_bar(stat = "identity",
           position = "dodge",
           color = "black") +
  facet_col(~ Z) +
  xlab("") +
  ylab("") +
  ggtitle("Spending on police services") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    plot.title = element_text(size = 14, face = "plain", hjust = 0.5),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 9),
    axis.title.x = element_text(size = 9),
    axis.text.y = element_text(size = 9),
    axis.title.y = element_text(size = 9)
  )


gg_df <-
  dat %>%
  select(Z,
         weights,
         othr_change_numcops_n) %>%
  rename(DV = othr_change_numcops_n) %>%
  mutate(
    DV = case_when(
      DV == 5 ~ "Increase a lot",
      DV == 4 ~ "Increase a little",
      DV == 3 ~ "Keep the same",
      DV == 2 ~ "Decrease a little",
      DV == 1 ~ "Decrease a lot"
    ),
    DV = factor(
      DV,
      levels = c(
        "Increase a lot",
        "Increase a little",
        "Keep the same",
        "Decrease a little",
        "Decrease a lot"
      )
    ),
    Z = factor(
      Z,
      levels = c("Abolish", "Defund", "Reform"),
      labels = c("Abolishing", "Defunding", "Reforming")
    )
  ) %>%
  filter(!is.na(DV)) %>%
  group_by(Z, DV) %>%
  summarise(n_wt = sum(weights),
            n = n()) %>%
  group_by(Z) %>%
  mutate(weighted = n_wt / sum(n_wt),
         unweighted = n / sum(n)) %>%
  pivot_longer(cols = weighted:unweighted,
               names_to = "type",
               values_to = "prop")

## Sanity check
gg_df %>% group_by(Z, type) %>% summarise(sum(prop))

g2 <-
  gg_df %>%
  ggplot(., aes(y = prop, x = DV, fill = type)) +
  geom_bar(stat = "identity",
           position = "dodge",
           color = "black") +
  facet_col(~ Z) +
  xlab("") +
  ylab("") +
  ggtitle("Number of police on the street") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    plot.title = element_text(size = 14, face = "plain", hjust = 0.5),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 9),
    axis.title.x = element_text(size = 9),
    axis.text.y = element_blank(),
    axis.title.y = element_text(size = 9)
  )


g1 <- g1 + theme(plot.margin = grid::unit(c(1, 0, 0, 0), "mm"))
g2 <- g2 + theme(plot.margin = grid::unit(c(1, 0, 0, 0), "mm"))


g <-
  plot_grid(
    g1,
    g2,
    labels = "auto",
    align = "h",
    label_fontface = "plain",
    label_fontfamily = "sans",
    label_size = 12
  )

g

g %>%
  ggsave(filename = "FigB8.1.pdf",
         .,
         width = 12,
         height = 7)

### ------ Fig. B8.2 Appendix ------

gg_df <-
  dat %>%
  select(Z_cut,
         weights,
         cuts_crime_n) %>%
  rename(DV = cuts_crime_n) %>%
  mutate(
    DV = case_when(
      DV == 5 ~ "Much higher",
      DV == 4 ~ "Slightly higher",
      DV == 3 ~ "About the same",
      DV == 2 ~ "Slightly lower",
      DV == 1 ~ "Much lower"
    ),
    DV = factor(DV,
                levels = rev(
                  c(
                    "Much higher",
                    "Slightly higher",
                    "About the same",
                    "Slightly lower",
                    "Much lower"
                  )
                )),
    Z_cut = factor(
      Z_cut,
      levels = c(1, 2, 3, 4),
      labels = c("10% cut", "25% cut", "50% cut", "100% cut")
    )
  ) %>%
  filter(!is.na(DV)) %>%
  group_by(Z_cut, DV) %>%
  summarise(n_wt = sum(weights),
            n = n()) %>%
  group_by(Z_cut) %>%
  mutate(weighted = n_wt / sum(n_wt),
         unweighted = n / sum(n)) %>%
  pivot_longer(cols = weighted:unweighted,
               names_to = "type",
               values_to = "prop")

## Sanity check
gg_df %>% group_by(Z_cut, type) %>% summarise(sum(prop))

g1 <-
  gg_df %>%
  ggplot(., aes(y = prop, x = DV, fill = type)) +
  geom_bar(stat = "identity",
           position = "dodge",
           color = "black") +
  facet_col(~ Z_cut) +
  xlab("") +
  ylab("") +
  ggtitle("Crime levels") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    plot.title = element_text(size = 14, face = "plain", hjust = 0.5),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 9),
    axis.title.x = element_text(size = 9),
    axis.text.y = element_text(size = 9),
    axis.title.y = element_text(size = 9)
  )


gg_df <-
  dat %>%
  select(Z_cut,
         weights,
         cuts_safe_n) %>%
  rename(DV = cuts_safe_n) %>%
  mutate(
    DV = case_when(
      DV == 1 ~ "Much higher",
      DV == 2 ~ "Slightly higher",
      DV == 3 ~ "About the same",
      DV == 4 ~ "Slightly lower",
      DV == 5 ~ "Much lower"
    ),
    DV = factor(DV,
                levels = rev(
                  c(
                    "Much higher",
                    "Slightly higher",
                    "About the same",
                    "Slightly lower",
                    "Much lower"
                  )
                )),
    Z_cut = factor(
      Z_cut,
      levels = c(1, 2, 3, 4),
      labels = c("10% cut", "25% cut", "50% cut", "100% cut")
    )
  ) %>%
  filter(!is.na(DV)) %>%
  group_by(Z_cut, DV) %>%
  summarise(n_wt = sum(weights),
            n = n()) %>%
  group_by(Z_cut) %>%
  mutate(weighted = n_wt / sum(n_wt),
         unweighted = n / sum(n)) %>%
  pivot_longer(cols = weighted:unweighted,
               names_to = "type",
               values_to = "prop")

## Sanity check
gg_df %>% group_by(Z_cut, type) %>% summarise(sum(prop))

g2 <-
  gg_df %>%
  ggplot(., aes(y = prop, x = DV, fill = type)) +
  geom_bar(stat = "identity",
           position = "dodge",
           color = "black") +
  facet_col(~ Z_cut) +
  xlab("") +
  ylab("") +
  ggtitle("Fear of crime") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_fill_manual("", values = c("white", "darkgrey")) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    plot.title = element_text(size = 14, face = "plain", hjust = 0.5),
    strip.background = element_blank(),
    legend.position = "none",
    axis.text.x = element_text(size = 9),
    axis.title.x = element_text(size = 9),
    axis.text.y = element_blank(),
    axis.title.y = element_text(size = 9)
  )


g1 <- g1 + theme(plot.margin = grid::unit(c(1, 0, 0, 0), "mm"))
g2 <- g2 + theme(plot.margin = grid::unit(c(1, 0, 0, 0), "mm"))


g <-
  plot_grid(
    g1,
    g2,
    labels = "auto",
    align = "h",
    label_fontface = "plain",
    label_fontfamily = "sans",
    label_size = 12
  )

g

g %>%
  ggsave(filename = "FigB8.2.pdf",
         .,
         width = 11,
         height = 8)

### ------ Sub-group estimates-----



### ------ Fig. B10a Appendix ----------
gg_u_df <-
  dat %>%
  select(supp_reform_n,
         supp_defund_n,
         supp_abolish_n,
         dem_race_4) %>%
  pivot_longer(supp_reform_n:supp_abolish_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(DV, dem_race_4) %>%
  do(tidy(lm_robust(Y ~ 1, data = .))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    DV = case_when(
      DV == "supp_reform_n" ~  "Reforming the police",
      DV == "supp_defund_n" ~  "Defunding the police",
      DV == "supp_abolish_n" ~ "Abolishing the police"
    ),
    type = "Slogan"
  ) %>%
  bind_rows(
    dat %>%
      select(
        policy_reform_n,
        policy_defund_n,
        policy_abolish_n,
        dem_race_4
      ) %>%
      pivot_longer(
        policy_reform_n:policy_abolish_n,
        names_to = "DV",
        values_to = "Y"
      ) %>%
      group_by(DV, dem_race_4) %>%
      do(tidy(lm_robust(Y ~ 1, data = .))) %>%
      filter(term == "(Intercept)") %>%
      mutate(
        DV = case_when(
          DV == "policy_reform_n" ~  "Reforming the police",
          DV == "policy_defund_n" ~  "Defunding the police",
          DV == "policy_abolish_n" ~ "Abolishing the police"
        ),
        type = "Substance"
      )
  )


gg_w_df <-
  dat %>%
  select(supp_reform_n,
         supp_defund_n,
         supp_abolish_n,
         dem_race_4,
         weights) %>%
  pivot_longer(supp_reform_n:supp_abolish_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(DV, dem_race_4) %>%
  do(tidy(lm_robust(Y ~ 1, weights = weights,  data = .))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    DV = case_when(
      DV == "supp_reform_n" ~  "Reforming the police",
      DV == "supp_defund_n" ~  "Defunding the police",
      DV == "supp_abolish_n" ~ "Abolishing the police"
    ),
    type = "Slogan"
  ) %>%
  bind_rows(
    dat %>%
      select(
        policy_reform_n,
        policy_defund_n,
        policy_abolish_n,
        dem_race_4,
        weights
      ) %>%
      pivot_longer(
        policy_reform_n:policy_abolish_n,
        names_to = "DV",
        values_to = "Y"
      ) %>%
      group_by(DV, dem_race_4) %>%
      do(tidy(lm_robust(
        Y ~ 1, data = ., weights = weights
      ))) %>%
      filter(term == "(Intercept)") %>%
      mutate(
        DV = case_when(
          DV == "policy_reform_n" ~  "Reforming the police",
          DV == "policy_defund_n" ~  "Defunding the police",
          DV == "policy_abolish_n" ~ "Abolishing the police"
        ),
        type = "Substance"
      )
  )

gg_df <-
  bind_rows(
    gg_u_df %>% mutate(estimator = "Unweighted"),
    gg_w_df %>% mutate(estimator = "Weighted")
  ) %>%
  mutate(type = factor(type, levels = c("Substance", "Slogan")))


tab_top <-
  gg_df %>%
  filter(estimator == "Weighted",
         dem_race_4 %in% c("Black", "White", "Hispanic")) %>%
  arrange(dem_race_4) %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -term, -outcome, -df, -estimator) %>%
  arrange(DV, desc(type), dem_race_4) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, type, dem_race_4, entry) %>%
  pivot_wider(names_from = type:dem_race_4,
              values_from = entry)

write.table(
  tab_top,
  file = "Tab10a_weighted.txt",
  sep = ",",
  quote = FALSE,
  row.names = F,
  col.names = FALSE
)

tab_top <-
  gg_df %>%
  filter(estimator == "Unweighted",
         dem_race_4 %in% c("Black", "White", "Hispanic")) %>%
  arrange(dem_race_4) %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -term, -outcome, -df, -estimator) %>%
  arrange(DV, desc(type), dem_race_4) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, type, dem_race_4, entry) %>%
  pivot_wider(names_from = type:dem_race_4,
              values_from = entry)

write.table(
  tab_top,
  file = "Tab10a_unweighted.txt",
  sep = ",",
  quote = FALSE,
  row.names = F,
  col.names = FALSE
)

g <-
  gg_df %>%
  filter(estimator == "Weighted", dem_race_4 != "Other") %>%
  ggplot(.,
         aes(
           x = estimate,
           y = type,
           group = dem_race_4,
           shape = dem_race_4,
           color = dem_race_4
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.7),
    size = 3,
    color = "white",
    fill = "black"
  ) +
  scale_shape_manual("Sub-group:", values = c(23, 22, 21)) +
  scale_color_manual("Sub-group:", values = c("black", "black", "black")) +
  facet_col( ~ DV) +
  scale_x_continuous(limits = c(1, 5)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "lines"),
    legend.margin = margin(
      t = -0.15,
      l = -0.30,
      unit = 'cm'
    ),
    legend.box.margin = margin(
      t = -0.15,
      l = -0.30,
      unit = "cm"
    ),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "bottom",
    legend.justification = "left",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Respondents' average level of support\n (1 = None at all,...,5 = A great deal)",
       y = "")

g <- g + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))

## Export results
g %>%
  ggsave(
    filename = "FigB10a.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * (nrow(g$data) * .5)
  )

### ------ Fig. B10b Manuscript ----------
gg_u_df <-
  dat %>%
  select(supp_reform_n,
         supp_defund_n,
         supp_abolish_n,
         dem_pid_3) %>%
  pivot_longer(supp_reform_n:supp_abolish_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(DV, dem_pid_3) %>%
  do(tidy(lm_robust(Y ~ 1, data = .))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    DV = case_when(
      DV == "supp_reform_n" ~  "Reforming the police",
      DV == "supp_defund_n" ~  "Defunding the police",
      DV == "supp_abolish_n" ~ "Abolishing the police"
    ),
    type = "Slogan"
  ) %>%
  bind_rows(
    dat %>%
      select(policy_reform_n,
             policy_defund_n,
             policy_abolish_n,
             dem_pid_3) %>%
      pivot_longer(
        policy_reform_n:policy_abolish_n,
        names_to = "DV",
        values_to = "Y"
      ) %>%
      group_by(DV, dem_pid_3) %>%
      do(tidy(lm_robust(Y ~ 1, data = .))) %>%
      filter(term == "(Intercept)") %>%
      mutate(
        DV = case_when(
          DV == "policy_reform_n" ~  "Reforming the police",
          DV == "policy_defund_n" ~  "Defunding the police",
          DV == "policy_abolish_n" ~ "Abolishing the police"
        ),
        type = "Substance"
      )
  )


gg_w_df <-
  dat %>%
  select(supp_reform_n,
         supp_defund_n,
         supp_abolish_n,
         dem_pid_3,
         weights) %>%
  pivot_longer(supp_reform_n:supp_abolish_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(DV, dem_pid_3) %>%
  do(tidy(lm_robust(Y ~ 1, weights = weights,  data = .))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    DV = case_when(
      DV == "supp_reform_n" ~  "Reforming the police",
      DV == "supp_defund_n" ~  "Defunding the police",
      DV == "supp_abolish_n" ~ "Abolishing the police"
    ),
    type = "Slogan"
  ) %>%
  bind_rows(
    dat %>%
      select(
        policy_reform_n,
        policy_defund_n,
        policy_abolish_n,
        dem_pid_3,
        weights
      ) %>%
      pivot_longer(
        policy_reform_n:policy_abolish_n,
        names_to = "DV",
        values_to = "Y"
      ) %>%
      group_by(DV, dem_pid_3) %>%
      do(tidy(lm_robust(
        Y ~ 1, data = ., weights = weights
      ))) %>%
      filter(term == "(Intercept)") %>%
      mutate(
        DV = case_when(
          DV == "policy_reform_n" ~  "Reforming the police",
          DV == "policy_defund_n" ~  "Defunding the police",
          DV == "policy_abolish_n" ~ "Abolishing the police"
        ),
        type = "Substance"
      )
  )

gg_df <-
  bind_rows(
    gg_u_df %>% mutate(estimator = "Unweighted"),
    gg_w_df %>% mutate(estimator = "Weighted")
  ) %>%
  mutate(type = factor(type, levels = c("Substance", "Slogan")))



tab_btm <-
  gg_df %>%
  filter(estimator == "Weighted") %>%
  arrange(dem_pid_3) %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -term, -outcome, -df, -estimator) %>%
  arrange(DV, desc(type), dem_pid_3) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, type, dem_pid_3, entry) %>%
  pivot_wider(names_from = type:dem_pid_3,
              values_from = entry)


write.table(
  tab_btm,
  file = "Tab10b_weighted.txt",
  sep = ",",
  quote = FALSE,
  row.names = F,
  col.names = FALSE
)

tab_btm <-
  gg_df %>%
  filter(estimator == "Unweighted") %>%
  arrange(dem_pid_3) %>%
  mutate_if(is.numeric, round, 3) %>%
  select(-statistic, -p.value, -term, -outcome, -df, -estimator) %>%
  arrange(DV, desc(type), dem_pid_3) %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error, digits = 2),
    DV = str_replace(DV, " the police", "")
  ) %>%
  select(DV, type, dem_pid_3, entry) %>%
  pivot_wider(names_from = type:dem_pid_3,
              values_from = entry)


write.table(
  tab_btm,
  file = "Tab10b_unweighted.txt",
  sep = ",",
  quote = FALSE,
  row.names = F,
  col.names = FALSE
)

g <-
  gg_df %>%
  filter(estimator == "Weighted") %>%
  ggplot(
    .,
    aes(
      x = estimate,
      y = type,
      group = dem_pid_3,
      fill = dem_pid_3,
      color = dem_pid_3,
      shape = dem_pid_3
    )
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(position = position_dodgev(height = 0.7),
             size = 3,
             color = "white") +
  scale_shape_manual("Sub-group:", values = c(23, 22, 21)) +
  scale_fill_manual("Sub-group:", values = bpr_colors) +
  scale_color_manual("Sub-group:", values = bpr_colors) +
  facet_col( ~ DV) +
  scale_x_continuous(limits = c(1, 5)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "lines"),
    legend.margin = margin(
      t = -0.15,
      l = -0.30,
      unit = 'cm'
    ),
    legend.box.margin = margin(
      t = -0.15,
      l = -0.30,
      unit = "cm"
    ),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "bottom",
    legend.justification = "left",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Respondents' average level of support\n (1 = None at all,...,5 = A great deal)",
       y = "")

g <- g + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g

## Export results
g %>%
  ggsave(
    filename = "FigB10b.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * (nrow(g$data) * .5)
  )


#### ----- Fig. B11a Manuscript -----
gg_u_df <-
  dat %>%
  select(Z,
         weights,
         othr_spend_police_n,
         othr_change_numcops_n,
         dem_race_4) %>%
  pivot_longer(
    othr_spend_police_n:othr_change_numcops_n,
    names_to = "DV",
    values_to = "Y"
  ) %>%
  group_by(Z, DV, dem_race_4) %>%
  do(tidy(lm_robust(Y ~ 1, data = .))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    DV = case_when(
      DV == "othr_spend_police_n" ~   "Spending on police services",
      DV == "othr_change_numcops_n" ~ "Number of police on the street"
    ),
    Z = case_when(
      Z == "Abolish" ~ "Abolishing",
      Z == "Defund" ~ "Defunding",
      Z == "Reform" ~ "Reforming"
    ),
    Z = factor(Z, levels = rev(
      c("Abolishing", "Defunding", "Reforming")
    ))
  )

gg_w_df <-
  dat %>%
  select(Z,
         weights,
         othr_spend_police_n,
         othr_change_numcops_n,
         dem_race_4) %>%
  pivot_longer(
    othr_spend_police_n:othr_change_numcops_n,
    names_to = "DV",
    values_to = "Y"
  ) %>%
  group_by(Z, DV, dem_race_4) %>%
  do(tidy(lm_robust(Y ~ 1, data = ., weights = weights))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    DV = case_when(
      DV == "othr_spend_police_n" ~ "Spending on police services",
      DV == "othr_change_numcops_n" ~ "Number of police on the street"
    ),
    Z = case_when(
      Z == "Abolish" ~ "Abolishing",
      Z == "Defund" ~ "Defunding",
      Z == "Reform" ~ "Reforming"
    ),
    Z = factor(Z, levels = rev(
      c("Abolishing", "Defunding", "Reforming")
    ))
  )

gg_top <-
  bind_rows(
    gg_u_df %>% mutate(estimator = "Unweighted"),
    gg_w_df %>% mutate(estimator = "Weighted")
  )

gg_u_df <-
  dat %>%
  select(Z_cut,
         weights,
         dem_pid_7n,
         cuts_crime_n,
         cuts_safe_n,
         dem_race_4) %>%
  pivot_longer(cuts_crime_n:cuts_safe_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(Z_cut, DV, dem_race_4) %>%
  do(tidy(lm_robust(Y ~ 1, data = .))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    estimator = "Unweighted",
    Z_cut = factor(Z_cut,
                   levels = rev(c(1, 2, 3, 4)),
                   labels = rev(c(
                     "10%", "25%", "50%", "100%"
                   ))),
    DV = case_when(
      DV == "cuts_crime_n" ~ "Crime levels",
      DV == "cuts_safe_n" ~ "Fear of crime"
    )
  ) %>%
  rename(Z = Z_cut)

gg_w_df <-
  dat %>%
  select(Z_cut,
         weights,
         dem_pid_7n,
         cuts_crime_n,
         cuts_safe_n,
         dem_race_4) %>%
  pivot_longer(cuts_crime_n:cuts_safe_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(Z_cut, DV, dem_race_4) %>%
  do(tidy(lm_robust(Y ~ 1, data = ., weights = weights))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    estimator = "Weighted",
    Z_cut = factor(Z_cut,
                   levels = rev(c(1, 2, 3, 4)),
                   labels = rev(c(
                     "10%", "25%", "50%", "100%"
                   ))),
    DV = case_when(
      DV == "cuts_crime_n" ~ "Crime levels",
      DV == "cuts_safe_n" ~ "Fear of crime"
    )
  ) %>%
  rename(Z = Z_cut)


gg_btm <- bind_rows(gg_u_df, gg_w_df)

gg_scale_df <-
  bind_rows(gg_top, gg_btm)  %>%
  select(-term)

tab_top <-
  gg_scale_df %>%
  filter(
    estimator == "Weighted",
    dem_race_4 != "Other",
    DV %in% c("Number of police on the street",
              "Spending on police services")
  ) %>%
  arrange(dem_race_4) %>%
  mutate_if(is.numeric, round, 3) %>%
  select(Z, DV, dem_race_4, estimate, std.error) %>%
  arrange(DV, dem_race_4) %>%
  mutate(entry = table_entry(est = estimate, se = std.error, digits = 2)) %>%
  select(Z, DV, dem_race_4, entry) %>%
  pivot_wider(names_from = DV:dem_race_4,
              values_from = entry)

write.table(
  tab_top,
  file = "Tab11a_weighted.txt",
  sep = ",",
  quote = FALSE,
  row.names = F,
  col.names = FALSE
)

tab_top <-
  gg_scale_df %>%
  filter(
    estimator == "Unweighted",
    dem_race_4 != "Other",
    DV %in% c("Number of police on the street",
              "Spending on police services")
  ) %>%
  arrange(dem_race_4) %>%
  mutate_if(is.numeric, round, 3) %>%
  select(Z, DV, dem_race_4, estimate, std.error) %>%
  arrange(DV, dem_race_4) %>%
  mutate(entry = table_entry(est = estimate, se = std.error, digits = 2)) %>%
  select(Z, DV, dem_race_4, entry) %>%
  pivot_wider(names_from = DV:dem_race_4,
              values_from = entry)

write.table(
  tab_top,
  file = "Tab11a_unweighted.txt",
  sep = ",",
  quote = FALSE,
  row.names = F,
  col.names = FALSE
)

tab_top <-
  gg_scale_df %>%
  filter(estimator == "Weighted",
         dem_race_4 != "Other",
         DV %in% c("Crime levels", "Fear of crime")) %>%
  arrange(dem_race_4) %>%
  mutate_if(is.numeric, round, 3) %>%
  select(Z, DV, dem_race_4, estimate, std.error) %>%
  arrange(DV, dem_race_4) %>%
  mutate(entry = table_entry(est = estimate, se = std.error, digits = 2)) %>%
  select(Z, DV, dem_race_4, entry) %>%
  pivot_wider(names_from = DV:dem_race_4,
              values_from = entry)

write.table(
  tab_top,
  file = "Tab12a_weighted.txt",
  sep = ",",
  quote = FALSE,
  row.names = F,
  col.names = FALSE
)

tab_top <-
  gg_scale_df %>%
  filter(
    estimator == "Unweighted",
    dem_race_4 != "Other",
    DV %in% c("Crime levels", "Fear of crime")
  ) %>%
  arrange(dem_race_4) %>%
  mutate_if(is.numeric, round, 3) %>%
  select(Z, DV, dem_race_4, estimate, std.error) %>%
  arrange(DV, dem_race_4) %>%
  mutate(entry = table_entry(est = estimate, se = std.error, digits = 2)) %>%
  select(Z, DV, dem_race_4, entry) %>%
  pivot_wider(names_from = DV:dem_race_4,
              values_from = entry)

write.table(
  tab_top,
  file = "Tab12a_unweighted.txt",
  sep = ",",
  quote = FALSE,
  row.names = F,
  col.names = FALSE
)

## Generate plots
g1 <-
  gg_scale_df %>%
  filter(str_detect(DV, regex("police", ignore_case = TRUE)),
         estimator == "Weighted",
         dem_race_4 != "Other") %>%
  ggplot(.,
         aes(
           x = estimate,
           y = Z,
           group = dem_race_4,
           shape = dem_race_4,
           color = dem_race_4
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.7),
    size = 3,
    color = "white",
    fill = "black"
  ) +
  scale_shape_manual("Sub-group:", values = c(23, 22, 21)) +
  scale_color_manual("Sub-group:", values = c("black", "black", "black")) +
  facet_col( ~ DV, space = "free") +
  scale_x_continuous(limits = c(1, 5)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "lines"),
    legend.margin = margin(
      t = -0.15,
      l = -0.30,
      unit = 'cm'
    ),
    legend.box.margin = margin(
      t = -0.15,
      l = -0.30,
      unit = "cm"
    ),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    legend.justification = "left",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Perceptions of supporters' preferences\n(1 = Increase a lot,...,5 = Decrease a lot)",
       y = "")

g2 <-
  gg_scale_df %>%
  filter(!str_detect(DV, regex("police|spending", ignore_case = TRUE)),
         estimator == "Weighted",
         dem_race_4 != "Other") %>%
  ggplot(.,
         aes(
           x = estimate,
           y = Z,
           group = dem_race_4,
           shape = dem_race_4,
           color = dem_race_4
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.7),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.7),
    size = 3,
    color = "white",
    fill = "black"
  ) +
  scale_shape_manual("Sub-group:", values = c(23, 22, 21)) +
  scale_color_manual("Sub-group:", values = c("black", "black", "black")) +
  facet_col( ~ DV, space = "free") +
  scale_x_continuous(limits = c(1, 5)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "lines"),
    legend.margin = margin(
      t = -0.15,
      l = -0.30,
      unit = 'cm'
    ),
    legend.box.margin = margin(
      t = -0.15,
      l = -0.30,
      unit = "cm"
    ),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    legend.justification = "left",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Perceived impact of X% cut to number of police\n(1 = Much lower,...,5 = Much higher)",
       y = "")

legend_b <- get_legend(
  g2 +
    guides(color = guide_legend(nrow = 1)) +
    theme(
      legend.position = "bottom",
      legend.margin = margin(t = -0.15, l = 0.30, unit = 'cm'),
      legend.box.margin = margin(t = -0.15, l = 0.30, unit = "cm"),
      legend.text = element_text(size = 10),
      legend.title = element_text(size = 10),
      legend.background = element_blank(),
      legend.box.background = element_blank(),
      legend.key = element_blank(),
      legend.justification = "left"
    )
)

g1 <- g1 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g2 <- g2 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))


pgrid <- plot_grid(
  g1,
  g2,
  ncol = 2,
  labels = "auto",
  label_fontface = "plain",
  label_fontfamily = "sans",
  label_size = 12
)

g <- plot_grid(pgrid,
               legend_b,
               ncol = 1,
               rel_heights = c(1, .1))
g

g %>%
  ggsave(filename = "FigB11a.pdf",
         .,
         width = 7.5,
         height = 6)


#### ----- Fig. B11b Manuscript -----
gg_u_df <-
  dat %>%
  select(Z,
         weights,
         othr_spend_police_n,
         othr_change_numcops_n,
         dem_pid_3) %>%
  pivot_longer(
    othr_spend_police_n:othr_change_numcops_n,
    names_to = "DV",
    values_to = "Y"
  ) %>%
  group_by(Z, DV, dem_pid_3) %>%
  do(tidy(lm_robust(Y ~ 1, data = .))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    DV = case_when(
      DV == "othr_spend_police_n" ~   "Spending on police services",
      DV == "othr_change_numcops_n" ~ "Number of police on the street"
    ),
    Z = case_when(
      Z == "Abolish" ~ "Abolishing",
      Z == "Defund" ~ "Defunding",
      Z == "Reform" ~ "Reforming"
    ),
    Z = factor(Z, levels = rev(
      c("Abolishing", "Defunding", "Reforming")
    ))
  )

gg_w_df <-
  dat %>%
  select(Z,
         weights,
         othr_spend_police_n,
         othr_change_numcops_n,
         dem_pid_3) %>%
  pivot_longer(
    othr_spend_police_n:othr_change_numcops_n,
    names_to = "DV",
    values_to = "Y"
  ) %>%
  group_by(Z, DV, dem_pid_3) %>%
  do(tidy(lm_robust(Y ~ 1, data = ., weights = weights))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    DV = case_when(
      DV == "othr_spend_police_n" ~ "Spending on police services",
      DV == "othr_change_numcops_n" ~ "Number of police on the street"
    ),
    Z = case_when(
      Z == "Abolish" ~ "Abolishing",
      Z == "Defund" ~ "Defunding",
      Z == "Reform" ~ "Reforming"
    ),
    Z = factor(Z, levels = rev(
      c("Abolishing", "Defunding", "Reforming")
    ))
  )

gg_top <-
  bind_rows(
    gg_u_df %>% mutate(estimator = "Unweighted"),
    gg_w_df %>% mutate(estimator = "Weighted")
  )

gg_u_df <-
  dat %>%
  select(Z_cut,
         weights,
         dem_pid_7n,
         cuts_crime_n,
         cuts_safe_n,
         dem_pid_3) %>%
  pivot_longer(cuts_crime_n:cuts_safe_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(Z_cut, DV, dem_pid_3) %>%
  do(tidy(lm_robust(Y ~ 1, data = .))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    estimator = "Unweighted",
    Z_cut = factor(Z_cut,
                   levels = rev(c(1, 2, 3, 4)),
                   labels = rev(c(
                     "10%", "25%", "50%", "100%"
                   ))),
    DV = case_when(
      DV == "cuts_crime_n" ~ "Crime levels",
      DV == "cuts_safe_n" ~ "Fear of crime"
    )
  ) %>%
  rename(Z = Z_cut)

gg_w_df <-
  dat %>%
  select(Z_cut,
         weights,
         dem_pid_7n,
         cuts_crime_n,
         cuts_safe_n,
         dem_pid_3) %>%
  pivot_longer(cuts_crime_n:cuts_safe_n,
               names_to = "DV",
               values_to = "Y") %>%
  group_by(Z_cut, DV, dem_pid_3) %>%
  do(tidy(lm_robust(Y ~ 1, data = ., weights = weights))) %>%
  filter(term == "(Intercept)") %>%
  mutate(
    estimator = "Weighted",
    Z_cut = factor(Z_cut,
                   levels = rev(c(1, 2, 3, 4)),
                   labels = rev(c(
                     "10%", "25%", "50%", "100%"
                   ))),
    DV = case_when(
      DV == "cuts_crime_n" ~ "Crime levels",
      DV == "cuts_safe_n" ~ "Fear of crime"
    )
  ) %>%
  rename(Z = Z_cut)


gg_btm <- bind_rows(gg_u_df, gg_w_df)

gg_scale_df <-
  bind_rows(gg_top, gg_btm)  %>%
  select(-term)

tab_btm <-
  gg_scale_df %>%
  filter(
    estimator == "Weighted",
    DV %in% c("Number of police on the street",
              "Spending on police services")
  ) %>%
  arrange(dem_pid_3) %>%
  mutate_if(is.numeric, round, 3) %>%
  select(Z, DV, dem_pid_3, estimate, std.error) %>%
  arrange(DV, dem_pid_3) %>%
  mutate(entry = table_entry(est = estimate, se = std.error, digits = 2)) %>%
  select(Z, DV, dem_pid_3, entry) %>%
  pivot_wider(names_from = DV:dem_pid_3,
              values_from = entry)


write.table(
  tab_btm,
  file = "Tab11b_weighted.txt",
  sep = ",",
  quote = FALSE,
  row.names = F,
  col.names = FALSE
)

tab_btm <-
  gg_scale_df %>%
  filter(
    estimator == "Unweighted",
    DV %in% c("Number of police on the street",
              "Spending on police services")
  ) %>%
  arrange(dem_pid_3) %>%
  mutate_if(is.numeric, round, 3) %>%
  select(Z, DV, dem_pid_3, estimate, std.error) %>%
  arrange(DV, dem_pid_3) %>%
  mutate(entry = table_entry(est = estimate, se = std.error, digits = 2)) %>%
  select(Z, DV, dem_pid_3, entry) %>%
  pivot_wider(names_from = DV:dem_pid_3,
              values_from = entry)


write.table(
  tab_btm,
  file = "Tab11b_unweighted.txt",
  sep = ",",
  quote = FALSE,
  row.names = F,
  col.names = FALSE
)

tab_btm <-
  gg_scale_df %>%
  filter(estimator == "Weighted",
         DV %in% c("Crime levels", "Fear of crime")) %>%
  arrange(dem_pid_3) %>%
  mutate_if(is.numeric, round, 3) %>%
  select(Z, DV, dem_pid_3, estimate, std.error) %>%
  arrange(DV, dem_pid_3) %>%
  mutate(entry = table_entry(est = estimate, se = std.error, digits = 2)) %>%
  select(Z, DV, dem_pid_3, entry) %>%
  pivot_wider(names_from = DV:dem_pid_3,
              values_from = entry)

write.table(
  tab_btm,
  file = "Tab12b_weighted.txt",
  sep = ",",
  quote = FALSE,
  row.names = F,
  col.names = FALSE
)

tab_btm <-
  gg_scale_df %>%
  filter(estimator == "Unweighted",
         DV %in% c("Crime levels", "Fear of crime")) %>%
  arrange(dem_pid_3) %>%
  mutate_if(is.numeric, round, 3) %>%
  select(Z, DV, dem_pid_3, estimate, std.error) %>%
  arrange(DV, dem_pid_3) %>%
  mutate(entry = table_entry(est = estimate, se = std.error, digits = 2)) %>%
  select(Z, DV, dem_pid_3, entry) %>%
  pivot_wider(names_from = DV:dem_pid_3,
              values_from = entry)

write.table(
  tab_btm,
  file = "Tab12b_unweighted.txt",
  sep = ",",
  quote = FALSE,
  row.names = F,
  col.names = FALSE
)

## Generate plots
g1 <-
  gg_scale_df %>%
  filter(str_detect(DV, regex("police", ignore_case = TRUE)),
         estimator == "Weighted",
         dem_pid_3 != "Other") %>%
  ggplot(.,
         aes(
           x = estimate,
           y = Z,
           group = dem_pid_3,
           fill = dem_pid_3,
           color = dem_pid_3
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.8),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.8),
    size = 3,
    color = "white",
    pch = 21
  ) +
  scale_fill_manual("Sub-group:", values = bpr_colors) +
  scale_color_manual("Sub-group:", values = bpr_colors) +
  facet_col( ~ DV, space = "free") +
  scale_x_continuous(limits = c(1, 5)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "lines"),
    legend.margin = margin(
      t = -0.15,
      l = -0.30,
      unit = 'cm'
    ),
    legend.box.margin = margin(
      t = -0.15,
      l = -0.30,
      unit = "cm"
    ),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    legend.justification = "left",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Perceptions of supporters' preferences\n(1 = Increase a lot,...,5 = Decrease a lot)",
       y = "")

g2 <-
  gg_scale_df %>%
  filter(!str_detect(DV, regex("police|spending", ignore_case = TRUE)),
         estimator == "Weighted",
         dem_pid_3 != "Other") %>%
  ggplot(.,
         aes(
           x = estimate,
           y = Z,
           group = dem_pid_3,
           fill = dem_pid_3,
           color = dem_pid_3
         )) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.8),
    size = 1.5
  ) +
  geom_point(
    position = position_dodgev(height = 0.8),
    size = 3,
    color = "white",
    pch = 21
  ) +
  scale_fill_manual("Sub-group:", values = bpr_colors) +
  scale_color_manual("Sub-group:", values = bpr_colors) +
  facet_col( ~ DV, space = "free") +
  scale_x_continuous(limits = c(1, 5)) +
  theme_fivethirtyeight(base_size = 14) +
  theme(
    plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "lines"),
    legend.margin = margin(
      t = -0.15,
      l = -0.30,
      unit = 'cm'
    ),
    legend.box.margin = margin(
      t = -0.15,
      l = -0.30,
      unit = "cm"
    ),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.text = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 10),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank(),
    strip.background = element_blank(),
    legend.position = "none",
    legend.justification = "left",
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(size = 10)
  ) +
  labs(subtitle = "",
       x = "Perceived impact of X% cut to number of police\n(1 = Much lower,...,5 = Much higher)",
       y = "")

legend_b <- get_legend(
  g2 +
    guides(color = guide_legend(nrow = 1)) +
    theme(
      legend.position = "bottom",
      legend.margin = margin(t = -0.15, l = 0.30, unit = 'cm'),
      legend.box.margin = margin(t = -0.15, l = 0.30, unit = "cm"),
      legend.text = element_text(size = 10),
      legend.title = element_text(size = 10),
      legend.background = element_blank(),
      legend.box.background = element_blank(),
      legend.key = element_blank(),
      legend.justification = "left"
    )
)

g1 <- g1 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))
g2 <- g2 + theme(plot.margin = grid::unit(c(-5, 1, 1, 1), "mm"))

pgrid <- plot_grid(
  g1,
  g2,
  ncol = 2,
  labels = "auto",
  label_fontface = "plain",
  label_fontfamily = "sans",
  label_size = 12
)

g <- plot_grid(pgrid,
               legend_b,
               ncol = 1,
               rel_heights = c(1, .1))
g

g %>%
  ggsave(filename = "FigB11b.pdf",
         .,
         width = 7.5,
         height = 6)